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c/3 To faithfully simulate ITER and other modern fusion devices, one must re- 
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R -1 solve electron and ion fluctuation scales in a five-dimensional phase space and time. 
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Simultaneously one must account for the interaction of this turbulence with the 
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slow evolution of the large-scale plasma profiles. Because of the enormous range of 

CNJ scales involved and the high dimensionality of the problem, resolved first-principles 
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global simulations are very challenging using conventional (brute force) techniques. 
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In this thesis, the problem of resolving turbulence is addressed by developing ve- 

O 

Q\ locity space resolution diagnostics and an adaptive collisionality that allow for the 

O 

^ confident simulation of velocity space dynamics using the approximate minimal nec- 

X 

essary dissipation. With regard to the wide range of scales, a new approach has been 
developed in which turbulence calculations from multiple gyrokinetic flux tube sim- 
ulations are coupled together using transport equations to obtain self-consistent, 
steady-state background profiles and corresponding turbulent fluxes and heating. 
This approach is embodied in a new code, Trinity, which is capable of evolv- 



ing equilibrium profiles for multiple species, including electromagnetic effects and 
realistic magnetic geometry, at a fraction of the cost of conventional global simula- 
tions. Furthermore, an advanced model physical collision operator for gyrokinetics 
has been derived and implemented, allowing for the study of collisional turbulent 
heating, which has not been extensively studied. To demonstrate the utility of the 
coupled flux tube approach, preliminary results from Trinity simulations of the 
core of an ITER plasma are presented. 
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temperature is significantly decreased by retaining kinetic electron 
effects 
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Chapter 1 
Introduction 



1.1 Motivation 

We are now approaching a significant milestone in the fusion program. Over 
the next eight years, a multi-billion dollar magnetic confinement device (the Inter- 
national Thermonuclear Experimental Reactor, or "ITER") will be built to demon- 
strate the feasibility of fusion as an alternative energy source. The design for this 
experiment reflects a myriad of advances made through experimental, theoretical, 
and numerical studies in our understanding of fundamental plasma processes. How- 
ever, there are still many issues critical to the success of ITER and to the economic 
and scientific feasibility of future fusion devices that are not well understood. 

The main goal of this thesis is to present a set of numerical tools and a sound 
numerical framework within which we can study one of the key fundamental physics 
issues for magnetic confinement fusion devices: the presence of anomalously high 
levels of particle, momentum, and energy transport observed in hot, magnetized 
plasmas. This anomalous transport, which is due to small-scale turbulence driven by 
localized instabilities (or "microinstabilities"), has been the subject of intense study 
within the fusion program for decades. Without this turbulence the performance of 
magnetic fusion devices would be considerably improved. For example, a turbulence- 
free Joint European Torus (JET) would reach fusion ignition. 

The presence of turbulence is certainly not inevitable. Indeed, JET, TFTR, 
DIII-D and other fusion devices have demonstrated operation with regions of the 
plasma essentially turbulence-free. Understanding, controlling, and ultimately re- 
ducing turbulence in magnetic fusion experiments is thus a formidable but achievable 
challenge for the fusion program. Much progress has been made in our qualitative 
understanding of turbulent transport, and in some cases quantitative agreement be- 
tween numerical simulations and experiment is remarkably good. However, plasma 
turbulence and equilibrium profile evolution are both complex problems, and first- 
principles simulations with experimentally relevant plasma parameters have by ne- 
cessity only addressed either the effect of turbulence on the equilibrium or vice 
versa. 

In this thesis, we present a rigorous theoretical and numerical framework that 
allows for the efficient simulation and routine study of the self-consistent interaction 
between plasma turbulence and equilibrium profiles. While our approach provides a 
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significant savings over direct global simulations, it is still very challenging numer- 
ically. We have therefore implemented velocity space resolution diagnostics and an 
adaptive collisionality that allow us to resolve simulations with an approximately 
minimal number of grid points in velocity space. Furthermore, our collision opera- 
tor is an improvement over previous operators: it possesses a number of desirable 
properties, including local conservation of particle number, momentum, and energy, 
and satisfaction of Boltzmann's if-Theorem. The latter property is of particular 
importance when considering equilibrium evolution, since it is necessary to ensure 
that system entropy is increased and equilibrium profiles are heated by collisions 
(instead of cooled). Our use of a theoretically sound collision operator also allows 
us to conduct quantitative studies of the effect of collisional heating on equilibrium 
profile evolution - a topic that has received little attention from the plasma physics 
community. 

We do not claim to have developed a numerical fusion device. There are 
a number of important processes currently neglected in our model (most notably 
the development of equilibrium shear flows and the physics of the edge pedestal, 
which critically affect the power output of fusion devices). However, the approach 
presented here provides a platform for studying novel effects that may arise from the 
self-consistent interaction between turbulence, transport, and heating. Furthermore, 
the code we have developed (named Trinity) is capable, within broad parameter 
ranges, of providing quantitative predictions of microstability thresholds, turbulent 
fluctuations and tokamak performance from first principles. 



The hot, magnetized plasmas present in magnetic confinement fusion devices 
are rich and complicated physical systems. They support an enormous spectrum 
of processes whose time and space scales span many orders of magnitude: heated 
to millions of degrees, charged particles spiral tightly around curved magnetic field 
lines at a significant fraction of the speed of light; the same particles drift slowly 
across magnetic field lines, transporting particles, momentum, and heat across the 
length of the device; a multitude of waves propagate through the plasma, from light 
waves to Alfven waves to drift waves; kinetic instabilities give rise to a sea of small- 
scale, rapidly fluctuating turbulence, and fluid instabilities can lead to bulk motion 
of the plasma and catastrophic disruptions. The time and space scales for some of 
the important processes that affect the performance of magnetic confinement fusion 



Each of these processes requires often complex modeling. Consequently, it is 
neither analytically nor numerically feasible to work with a single model that simul- 
taneously describes all of the physical processes present. Instead, we must determine 
which processes are of greatest interest and identify reasonable approximations that 
will allow us to develop simplified models of their behavior. Occasionally, we may 
gain insights from these simplified numerical models that allow further reductions 
of the problem, but this cannot always be achieved because of the large number of 



1.2 Multiple scales 
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parameters and interactions that are known to be important experimentally. 

There are many important issues that must be addressed in order to develop 
a scientifically and economically viable fusion reactor. Fundamentally, however, we 
are interested in achieving high core pressures with minimal power input. This 
requires minimizing the radial heat transport, which is due primarily to turbulence. 

In order to address the challenges associated with turbulent transport, it is 
generally believed that one must take into account the close coupling between the 
slow (~ 1 s) evolution of large-scale (~ 1 m) variations in equilibrium density, 
temperature, and flow profiles and the rapid (~ 1 MHz) fluctuations of small- 
scale (~ m) plasma turbulence. This interaction of vastly disparate temporal 
and spatial scales renders direct numerical and analytical approaches intractable; 
instead, more sophisticated multiscale models are required. One such model is 
derived in Chapter [3j with the notable absence of equations describing the evolution 
of equilbrium flows. These flows are believed to play a critical role in the formation 
of the edge pedestal and internal transport barriers, thus limiting the immediate 
applicability of Trinity to core plasmas. 

To overcome the difficulty associated with the presence of a wide range of 
scales, different models have typically been applied to address turbulence and trans- 
port separately [H El El lU El El E|- Slowly-evolving, large-scale plasma transport 
has been widely modeled as a diffusive process, with theoretically and numerically 
derived diffusion coefficients. The magnetic equilibrium is typically modeled with 
the equations of magnetohydrodynamics (MHD), which treat the plasma as a single 
magnetized fluid. These approaches are generally inadequate for accurately describ- 
ing the rapidly-evolving, small-scale turbulence responsible for anomalous transport 
in fusion devices. The instabilities driving microturbulence arise, in part, due to 
the development of nontrivial structure in the distribution of plasma particle veloc- 
ities (which can be present due to the long collisional mean free path in hot fusion 
plasmas). Since this structure is not easily captured by conventional fluid models 
or tractable analytical approaches, a numerical description of kinetic, small-scale 
plasma turbulence is necessary. An example illustrating this point is provided in 
Chapter |2} 

1.3 Kinetic nature of magnetized plasma tur- 
bulence 

In order to address the complexities of plasma turbulence with existing com- 
puter technology, the full kinetic description must be simplified. This can be ac- 
complished by exploiting the separation of time and space scales in fusion plasmas. 
In this thesis, we employ the widely-used 5f gyrokinetic model El ED], which 
takes advantage of the following scale separations: the turbulence and resultant 
fluxes are calculated in a stationary equilibrium, exploiting the separation of the 
fast turbulence time scale and the slow profile evolution time scale; the variation of 
equilibrium gradient scale lengths perpendicular to the magnetic field line is ignored 
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Physics 


Space scale 


Time Scale 


Electron Energy 
Transport from 
ETG modes 


Scale perpendicular to B is 
~ Pe — Pi ~ 0.001 cm — 0.1 cm 
Scale parallel to B is 
qR ~ 15 m 


uj* e ~ 500 kHz - 5 MHz 


Ion Energy 
Transport from 
ITG modes 


Scale perpendicular to B is 
~ Pi — V ' PiLr ~ 0.1 cm — 8 cm 
Scale parallel to B is qR ~ 15 
m 


u* ~ 10- 100 kHz 


Transport 
Barriers 


Unknown scaling of 
perpendicular scales. 
Measured scales suggest width 
~ 1 — 10 cm 


Lifetime 100 s or more in 
core? Relaxation oscilla- 
tions for edge barrier with 
unknown frequency. 


Magnetic 
islands, 

Tearing modes 
and NTMs. 


Island width ~ lOp^ ~ 1 cm. 

Eigenfunction extent ~ L p ~ 
100 cm. Turbulent correlation 
length near island ~ 1 cm? 


Growth time ~ 1 — 100 s. 

Island frequency ~ 100 Hz — 
1 kHz. Turbulent frequency 
near island ~ 100 kHz 



Table 1.1: Some important tokamak space and time scales. Numerical values refer 
to ITER. 



(local assumption), exploiting the separation of the short perpendicular turbulence 
scale and the long perpendicular profile scale; and the dynamics of the turbulence 
itself is calculated assuming the particles gyrate about the ambient magnetic field 
lines infinitely fast, exploiting the difference in time scales between the dynamics of 
interest and a host of much faster processes that occur in magnetized plasmas. Fur- 
thermore, a distinction is made between fluctuations along the equilibrium magnetic 
field, which are assumed to have long (device size) wavelengths, and cross-field fluc- 
tuations, which have short (Larmor radius) wavelengths. Finally, the experimentally 
observed and theoretically well-founded expectation that the turbulent correlation 
lengths in the directions perpendicular to the magnetic field are small compared to 
the device dimensions (for large enough devices, high enough magnetic fields, and 
suitable distances from edge boundary layers) allows one to simulate small volumes 
of plasma surrounding individual magnetic field lines, called flux tubes, and to ex- 
trapolate the results from these small volumes to nearby flux tubes [TT]. [See Sec. 



1.6 for a more detailed discussion.] This is an assumption of statistical homogeneity 



among patches of plasma that are many turbulent correlation lengths apart. It is 
a particularly well-motivated and unsurprising approach for axisymmetric confine- 
ment devices such as tokamaks. It would be unwise to ignore this opportunity to 
reduce the simulation effort, choosing instead to simulate a large number of sta- 
tistically identical regions of plasma, absent an expectation of something such as 
important intermittent fluctuations. 
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Physics 


Space scale 


Time Scale 


Sawteeth 


Reconnection layer width 
~ 0.05 cm 

Eigenf unction extent ~ L p ~ 
100 cm. 


Crash time 50 fis — 100 [is 
Real frequency ~ 100 Hz — 
lkhz. 

Ramp time 1 — 100s 


Discharge 
Evolution 


Profile scales L p ~ 100 cm 


Energy confinement time 2— 
4s 

Burn time unknown 



Table 1.2: Some important tokamak space and time scales. Numerical values refer 
to ITER. 




Figure 1.1: Illustration of the flux tube simulation domain used in Trinity. Colors 
represent the amplitude of perturbations in the electrostatic potential. Notice that 
the turbulence is long wavelength along the equilibrium magnetic field and short 
wavelength in the plane perpendicular to it. Graphic courtesy of D. Applegate. 



These assumptions allow for the reduction of the problem from the long-time 
evolution of fast, gyroradius-scale turbulence throughout the full device, to the slow 
evolution of a few coupled magnetic flux tubes, each filled with fast, small-scale 
turbulent fluctuations. The fundamental validity of this approach for sufficiently 
large device size (p* ~ 0.003) has been demonstrated [12J by comparing results 
from flux tube simulations with results for the same cases from global simulations 



(Fig. 1.2 ), which allow for radial variation of equilibrium profiles within a turbulence 



simulation. 

Despite the significant simplifications granted by these gyrokinetic assump- 
tions, plasma turbulence simulations are still computationally challenging. Tur- 
bulence in conventional, neutral fluids is already a complex phenomenon; under- 
standing it has proven to be one of the great scientific challenges of our time. Ki- 
netic plasma turbulence, which may be characterized as particles interacting pri- 
marily with electromagnetic waves and occasionally with one another via collisions, 
possesses an additional level of complexity. For instance, a fundamental concept 
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0.2 0.4 0.6 0.8 

r/a 

Figure 1.2: Comparison of ion thermal diffusivity Xi calculated from local (GS2) and 
global (GYRO) simulations as a function of p* = p/a, where p is the gyroradius and 
a is the minor radius of the device. For sufficiently small p*, the local and global 
calculations of thermal diffusivity are in excellent agreement. Figure taken from 
Ref. [12]. 
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Figure 1.3: (Left): Cascade of entropy from large to small physical space struc- 
tures (W* ~ E|k ± |=fc x ^o|$k| 2 /2T and W h ~ E\^\= k± I d ^ T \h k \ 2 /2F are 
the entropy generation arising from the Boltzmann and non-Boltzmann responses 
of the perturbed distribution function, respectively). Solid black lines are theoreti- 
cal predictions [15J, and colored lines are data taken from 4D (k\\ = 0), electrostatic 
turbulence simulations at different resolutions [16]. (Right): Spectra characteriz- 
ing the cascade of entropy in velocity space, with E g (p) = J2kP l#k(p)| 2 , where 
9k(p) — f d 3 v Jo(pv±)g] l (y) is the Hankel transform of the guiding center perturbed 
distribution function, g. Solid black line is the theoretical prediction [IB], and col- 
ored lines are data taken from same runs as the figure on the left. Figures taken 
from Ref. [IB"] . 



in fluid turbulence is the cascade of energy from large-scale to small-scale spatial 
structures. In gyrokinetic turbulence, the three-dimensional cascade is replaced by 
a five-dimensional cascade of entropy from large-scale to small-scale phase space 



structures [El EH E51 GE1 021 ■ This is illustrated in Fig. 1.3 



It is well known that in weakly collisional plasmas, Landau and Barnes damp- 
ing of electromagnetic fluctuations leads to the development of small-scale structure 
in the distribution of particle parallel velocities. This is a result of mixing in phase 
space, where particles streaming along magnetic field lines transfer spatial structure 
into velocity structure. [The potential development of infinitesimally small scales 
in velocity space is illustrated in Appendix [Bj where we consider the simple case of 
collisionless Landau damping of the ion acoustic wave.] In addition to this paral- 
lel phase mixing arising from linear convection, there exists a perpendicular phase 
mixing process due to the averaged E x B particle motion [13J. A cartoon of this 



process is shown in Fig. 1.4 As particles rapidly gyrate about equilibrium magnetic 
field lines, they see spatially varying electromagnetic fluctuations that are essentially 
static in time. The E x B drift they experience is thus a result of the gyroaver- 
aged electromagnetic fields they see. Particles at the same guiding center position 
experience different gyroaveraged fields depending on their Larmor radius (which 
depends on perpendicular particle velocities) and thus drift with different guiding 
center velocities. This results in a mixing of particles with different perpendicular 
velocities in the gyrocenter distribution function and the generation of small scales 
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in the perpendicular velocity space. 

Because of the tendency of weakly-collisional plasmas to develop fine struc- 
tures in velocity space, one must pay careful attention to numerical resolution in 
velocities. This is ideally done by conducting a grid convergence study in velocity 
space. However, this is a numerically expensive process, so it is not always done. 
We have developed computationally cheap velocity space resolution diagnostics that 
allow us to monitor resolution at virtually no additional cost. When coupled with an 
adaptive collisionality, we are able to confidently simulate velocity space dynamics 
with the approximate minimal number of grid points in velocity space necessary for 
resolution. This is detailed in Chapter 4. 

1.4 Turbulent heating and the importance of 
collisions 

The evolution of equilibrium pressure profiles is determined by a balance be- 
tween transport processes and local heating. The net heating consists of contribu- 
tions from a number of sources, including external heating, atomic heating, Ohmic 
heating, and thermal energy exchange between species. Most of these phenomena 
have been extensively studied, both analytically and through the use of numerical 
transport solvers. However, little attention has been given to anomalous heating 
arising from microturbulence. The gyrokinetic turbulent heating of species s can be 
defined in various ways. In Refs. [IH] and [20] it is defined to be 



where 83 \\ is the perturbed parallel current, 83d is the current perturbation due 
to particle drifts, and 5Eii and <5Ej_ are the parallel and perpendicular perturbed 
electric fields, respectively. In Chapter 3, we derive an equation for the evolution 
of the equilibrium pressure that leads us to a somewhat different definition for the 
turbulent heating. We show, however, that both definitions lead to a net (species- 
summed) turbulent heating of zero. 

While the net turbulent heating is zero, the turbulent heating for each species 
(or equivalently, the turbulent energy exchange between species) is not necessarily 
zero. It is formally the same order as the heat transport, so there is a possibility 
of significant turbulent energy exchange between species. In the cases considered 
in Ref. [20J, it was found that the parallel and perpendicular contributions to the 
turbulent heating nearly cancel, giving only a 10% adjustment to the net heating. 
However, to our knowledge, no additional cases have been considered, and the def- 
inition used for the turbulent heating does not contain all of the turbulent heating 
terms appearing in the equations we derive in Chapter 3 for the time evolution of 
the equilibrium pressure. Turbulent heating therefore deserves further careful study. 

In Chapter 3, we also express the turbulent heating as the sum of a positive- 
definite quantity describing collisional entropy generation and a term representing 
energy exchange between the equilibrium and the turbulence. Because the collisional 
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Figure 1.4: Cartoon illustrating nonlinear perpendicular phase mixing, which leads 
to the development of small-scale structures in both physical and velocity space. 
When the separation between particle gyroradii with the same guiding center be- 
comes comparable to the characteristic wavelength of the turbulence, the motion of 
the particles become decorrelated. Since the size of the gyroradii are proportional 
to the particles' perpendicular velocities, this indicates a decorrelation in velocity 
space as well, leading to the development of small-scale structure. Figure taken from 
Ref. [IE]. 
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entropy generation term is positive definite, it is generally easier to obtain a con- 
verged statistical average for it in numerical simulations than for the 53 ■ terms, 
which tend to have large amplitude oscillations associated with particles "sloshing" 
back and forth in plasma waves. To calculate the collisional entropy generation, we 
have developed a model gyrokinetic collision operator that retains the key properties 
of physical collisions. 

In general, collisional physics are not carefully treated (if treated at all) in 
gyrokinetic simulations of turbulence [2"Tj 122] . In principle, one would like to in- 
clude the full linearized Landau operator [22], but numerical implementation and 
calculation of the so-called "field-particle" part of this operator is quite challenging. 
Approximate models for the collision operator have been derived [2H [25] . but they 
do not possess all of the properties one would like a collision operator to possess [21] . 
Consequently, we have derived a new model collision operator for use in gyrokinet- 
ics, which is an improvement over previous operators. Some of the key properties 
of our operator are: local conservation of particle number, momentum, and energy; 
satisfaction of Boltzmann's if-Theorem; efficient smoothing in velocity space; and 
reduction to the full linearized Landau operator in the short wavelength limit (where 
dissipation primarily occurs). The derivation of this operator is presented in Chap- 
ter 5, and numerical implementation in Trinity and tests are presented in Chapter 
6. 



1.5 Stiff transport 

Kinetic microinstabilities depend on a large number of plasma parameters. 
However, the dominant microinstabilities in most magnetic confinement fusion de- 
vices are driven unstable primarily by sufficiently strong temperature gradients. 
Since these microinstabilites cause high levels of heat transport, they effectively 
limit temperature gradients to be at or below the critical gradient at which the 
kinetic modes go unstable (unless the temperature near the edge of the plasma is 
low or the external heating is very large) [261 [271 [28] . Below the critical gradi- 
ent, there is a low level of transport due to neoclassical effects (see e.g. Ref. [29J) 
that has a relatively weak dependence on temperature gradient scale length. Above 
the critical gradient, the level of transport increases dramatically because turbulent 
transport has a stiff dependence on temperature gradient scale length. For relatively 
high temperature plasmas with reasonable external heating power, this feature (stiff 
transport) leads to profiles adjusting so that their gradients are stuck at the critical 
gradient. Consequently, the core temperature depends sensitively on the tempera- 



ture at the edge of the device (Fig. [L5). Without high edge temperatures, ITER will 
not likely achieve its target core temperature, for instance [30]. The edge plasma is 
not modeled in this thesis because of the complicated physics involved and because 
sharp gradients occur near the edge of the device (in what is known as the edge 
pedestal), challenging the applicability of the gyrokinetic ordering we consider. 

This stiff dependence of the fluxes on the driving gradients and the sharp 
transition between neoclassical and turbulent transport at the critical gradient have 
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TW ( P| -0.8) [keV] 

Figure 1.5: Plot of the core temperature as a function of the edge temperature in 
the High-confinement mode of operation (H-mode) on ASDEX-U. Note the linear 
scaling, which indicates that the temperature gradient scale length across the device 
is fixed (at the critical gradient) and independent of temperature. The tendency of 
profile gradients to stay near the critical gradient implies a stiff dependence of the 
heat flux on equilibrium gradients. Figure taken from Ref. [HTj . 



11 



another unfortunate consequence: they make turbulent transport simulations very 
challenging. Stiff systems are notoriously difficult to address numerically because 
the sensitivity of the equations to small perturbations can lead to extreme restric- 
tions on the time step size. In order to avoid (or at least limit) these restrictions, one 
should treat the transport equations implicitly. Developing such an implicit scheme 
is a nontrivial problem since the transport is described by a set of coupled, nonlinear 
partial differential equations. However, implicit techniques for nonlinear equations, 
such as Newton's method, have successfully been applied to plasma transport equa- 
tions with model fluxes [52] ■ We derive an implicit technique for solving the plasma 
transport equations with nonlinear, gyrokinetic fluxes in Chapter 7. 



1.6 Multiscale simulations of turbulent trans- 
port and heating 

Assuming no intermediate time or space scales are present, a direct numerical 
simulation resolving fine (turbulence) time and space scales throughout the volume 
of a fusion device for an entire discharge is not necessary. Instead, one can use 
the separation of scales embodied in the gyrokinetic turbulence and transport equa- 
tions derived in Chapter [3] to embed small regions of fine grid spacing in a coarse, 



equilibrium-scale mesh (Fig. 1.6). We adopt this approach by calculating turbulent 



fluxes and heating in a series of flux tubes, each of which is used to map out an 



entire magnetic flux surface (Fig. 1.7). These flux surfaces are coupled together as 
radial grid points in the one- dimensional equations describing the evolution of radial 
profiles of equilibrium density and pressure. 

The computational savings from using our multiscale scheme can be quite 
large. The use of field line-following coordinates decreases the number of grid points 
necessary along the equilibrium magnetic field since parallel turbulence wavelengths 
are much longer than perpendicular wavelengths. The use of a flux tube simulation 
domain to map out an entire flux surface decreases the number of grid points nec- 
essary in the direction perpendicular to the field (but lying near the flux surface). 
Although the radial domain covered by a series of coupled flux tubes is comparable 
to the domain of a conventional global approach, the spacing of the radial grid points 
is more optimal. This is because the range of wavenumbers (or equivalently, the grid 
spacing) necessary to resolve the turbulent fluctuations varies across the large-scale 
radial profile due to variations in density, temperature, and magnetic geometry. 
Each flux tube is naturally able to simulate a range of wavenumbers independent of 
the other flux tubes, constituting an adaptive radial grid. Finally, evolution of the 
turbulence and transport on separate time scales using the gyrokinetic hierarchy of 
Chapter 3 allows for simulation of the entire discharge while sampling only a frac- 
tion of the total discharge time. [Note that the algorithms derived and implemented 
here can be used to simulate the time-dependent evolution of the equilibrium, even 
for "fast" phenomena, such as heat and cold pulses; steady-state transport is not 
assumed.] Taking into account all of these contributions, the rough savings estimate 
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Figure 1.6: (Center): Fine scale grid in space and time. (Top left): Coarse equi- 
librium grid spacing in time, with regions of fine grid spacing embedded. Each 
horizontal red strip represents simulation of turbulent dynamics to steady-state, 
keeping equilibrium quantities constant. (Top right): Coarse equilibrium grid spac- 
ing in radius, with regions of fine grid spacing embedded. Each vertical green strip 
represents simulation of turbulent dynamics in a narrow flux tube, assuming no ra- 
dial variation of equilibrium profiles or gradients across the domain. (Bottom left): 
Combination of the multiscale space and time grids. (Bottom right): Small blue 
squares are the simulation domain resulting from the multiscale mesh in space and 
time. 
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Figure 1.7: Illustration of flux tubes from Trinity simulations. Using statistical 
periodicity of the turbulence, a single flux tube (top left) several decorellation lengths 
long can be used to map an entire flux surface (3 flux tubes at top right, 6 at bottom 
left, and 8 at bottom right). Colors represent the amplitude of perturbations in the 
electrostatic potential. Graphics courtesy of D. Applegate. 
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given for ITER in Chapter 7 is a factor on the order of 10 10 . These savings can be 
used to include additional physics, such as coupled electron-ion dynamics, electro- 
magnetic fluctuations, multiple ion species, etc., in each flux tube. Furthermore, 
they place coupled turbulence, transport, and heating calculations within reach on 
current computing resources. 

At the time of this writing, it is possible to obtain millions of CPU-hours on 
parallel computers with O(10 5 ) processors. For a global simulation of a steady-state 
ITER core plasma, Trinity might require 16 flux tubes, each running turbulence 
simulations requiring ~ 4000 processors - enabling multispecies, electromagnetic 
turbulence simulations in each flux tube, for example. The algorithm derived below 
can spawn 2-4 copies of each flux tube simultaneously to estimate the fluxes and their 
main dependencies; the precise number can be determined at run time to match the 
available resources. Assuming 2 copies of each of the 16 flux tubes, each running on 
4000 processors, such a simulation would utilize 128,000 cores, with nearly perfect 
linear scaling, and should run to completion in a few hours. Thus, the algorithms 
presented here will allow routine simulations to study a range of physical conditions 
and magnetic configurations on existing computers, not just an annual "stunt run" 
with limited physics content and limited scientific value. 
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Chapter 2 
Microstability 



2.1 Introduction 

Kinetic theory is complicated, but sometimes necessary. In the hot, magne- 
tized plasmas of magnetic confinement fusion experiments, the collisional mean free 
path can be many kilometers - distances much greater than the device size. This 
leads to the development of nontrivial structure in the distribution of particle ve- 
locities, as we discuss in detail in Chapter |4j Conventional fluid models do not 
accurately describe drift-type instabilities under these circumstances [33]. Because 
drift instabilities typically induce strong energy transport when the driving gradient 
is pushed beyond the threshold of the given instability, knowledge of the threshold 
criterion is key to the interpretation of much experimental data. 

In this chapter we illustrate the necessity of a kinetic treatment for the insta- 
bilities leading to small-scale plasma turbulence. We do so by calculating a kinetic 
stability threshold and comparing with stability thresholds from various fluid theo- 
ries [3U |35l |36l [37]. What we will find is that fluid theory significantly underesti- 
mates the range of instability [38J. 

As our example system, we choose to consider the entropy mode [39J in a 
Z-pinch magnetic field configuration [40J. This configuration consists of a current 
running through the plasma in the z direction, generating a radially varying equilib- 
rium magnetic field in the </> direction. Here, we are using cyclindrical coordinates, 
i.e. (R,<f>,z). For strong pressure gradients, the plasma is unstable to magneto- 
hydrodynamic (MHD) instabilities with fast growth rates [Ml SI]- If the pressure 
gradient is sufficiently weak, the plasma is stable to MHD instabilities, but poten- 
tially unstable to the entropy mode. To demonstrate the importance of the kinetic 
approach, we calculate the stability threshold of the low-/? (electrostatic) entropy 
mode and compare with the results obtained from a number of fluid theories. 

2.2 Linear stability analysis 

For our linear stability analysis, we will be working within the framework of 
5f gyrokinetics, which is described in detail in Chapter [3] The distribution function 
/ for species s is given by 

f s = F 0s + h a + f 2a , (2.1) 
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where h s is the non-Boltamann part of the lowest order perturbed distribution func- 
tion, f 2s contains higher order terms, and F 0s = F Ms (1 — q s $/T 0s ), with F M a 
Maxwellian, q s the particle charge, $ the electrostatic potential, and T 0s the equi- 
librium temperature. With these definitions, the electrostatic version of the linear, 
collisionless gyrokinetic equation is 

^ + «„b ■ Vh s + (y E ) H ■ VF 0s + v B • Vh s = g ;^($)R (2 2) 

where b = B /-Bo is the unit vector in the direction of the equilibrium magnetic 
field, B , 

w E = — b x V$ (2.3) 
Bo 



is the E x B velocity, 



v B = — x [vf, + — - - -5-^? J_l (2.4) 



fi V " 2 ; 5 c5 fi ■■ 

is the sum of the curvature and VB drift velocities, f2o = qBo/mc is the particle 
gyrofrequency, is the equilibrium perpendicular current, and the angled brackets 
(.} R denote an average over gyroangle at fixed guiding center position R. 

To proceed, we use the form of the equilibrium magnetic field, Bo = Bo(r)<p, to 
compute 3± and to determine an MHD equilibrium condition. After some algebra, 
we find 

4 V = ^(l + #S) (2-5) 



c R V. B OR 

VBo 1 A „ R 



B HK 1 -^)*' < 2 ' 6 ' 

where (3 = 8ttp /Bq is the plasma beta, p is the equilibrium pressure, and L~ x = 
—dlnpo/dr is the inverse pressure gradient scale length. In the low (3 limit, we find 
dBo/dR ps —B /R, giving Jj_ ps and 



We now return to the linear gyrokinetic equation (2.2). For simplicity, we 
take the ion and electron temperature gradients to be zero. Assuming perturbed 
quantities are of the form h = hexp[ik z z — iut), we obtain an algebraic equation for 
h: 



k, ( o v 



Rn 0s V 11 2 



h s = - (u - lu* s ) F 0s , (2. 

-LOs 



where 



ck z T 0a 

u * s = a n t (2 ' 9) 
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is the diamagnetic drift frequency. Defining the normalized quantities u>n s = u/^*s 
and x = v/vth, a an d solving for h s , we have 



T o* ^ - ^n[g s ] |X n /i2| (xj s + x 2 ±s /2 ) T 0s /T ( 



(2.10) 



' Oe 



We currently have an additional unkown: $. We can obtain an expression for $ by 
using Poisson's equation and asserting quasineutrality (i.e. Y^ s 1s n s = 0). In terms 
of h and quasineutrality gives 



n.s 



0. 



where (.) r represents a gyroaverage at constant particle position r. 



Substituting Eq. (2.10) into Eq. (2.11), we obtain 



(2.11) 



E 



d 3 v J (a s 



Wat s - Sgn[q s ]T Qs /T ( 



Oe 



lu Ns - Sgn[q s ] \L n /R\ (x 2 + x 2 ±s /2) T 0s /T ( 



1 F 



Ms 



Oe 



(2.12) 

where Jq is a Bessel function of the first kind, and a = k z v±/Qo. The velocity 



integration in Eq. (2.12) is nontrivial. To simplify the analysis, we focus on the 



limit in which k±p s <C 1 (the drift-kinetic limit). In this case, Jo(a) « 1. We are 
then interested in evaluating an integral of the form 



exp[— x 2 



u N - Sgn[q s }£ ( xjj + x\j2 



(2.13) 



where £ = \L n /R\(T 0s /T 0e ). Defining Co = u N /£+Sgn[q s ]x 2 ± /2 and I N = -I£/2nv^ s , 
and using cylindrical velocity space coordinates, our integral takes the form 



'TV 



dx± x± exp[— x 2 



dx 



'o J-oc "w + Sgn[q a ]x^ 

Focusing on the x» integral, we take an aside and consider 

exp[— asj] 

dxw 



exp[— x 2 ] 



2 • 



(2.14) 



io + Sgn[q s ]x 



2 • 



(2.15) 



where a is a velocity-independent parameter. Differentiating I with respect to a 
gives 

dl f°° , 2 exp[-axj] 

— i CajQu II JC 



(2.16) 
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Now we integrate this expression by parts: 
dl 



da 



Sgn[q s 



Sgn[q s ] QI 



dx\\ exp[— ocxii] + Co I dx\\ — 



exp[— ax\ 



Co + Sgn[q s ]xf, 



(2.17) 



This differential equation has two distinct solutions depending on the value of 
Sgn[q s ). They are 



J_ = exp[— aw] I C\ + 7r 



I + = exp[au;] I C 2 — 7r 




(2.18) 



(2.19) 



where the subscript on I denotes the sign of q s , Erf is the error function, Erfi is the 
imaginary error function, and C\ and Ci are unspecified constants. 



In order to determine C\ and C*2, we must go back to Eq. (2.15), set a = 0, 
and perform the resulting integral. We find 



I{a = 0) 



da; 1 1 (uj + Sgn[q s ]xT) 1 



7T 



Sgn[q, 



10 



(2.20) 



Once again, this splits into two solutions depending on the value of Sgn[q s ]: 



I_(q! = 0) = Sgn[Im[uN]]i 



71 



10 



I+{a = 0) 



7T 



— ■ 



(2.21) 
(2.22) 



where Im[UN] denotes the imaginary part of a; at. Applying these results to Eqs. (|2.18|) 

>wing expressions fc 

Ci = Sgn[Im\uj]s[\\i 



and (2.19), we obtain the following expressions for G\ and 

7T 



Co 



71 



The solutions for /_ and I + are then 
J_ = exp[— au} 



Sgn\Im[ujfq\}i + Erfifv^CM^I 



71 



J + = expfacj]— =Erfc[v ad)] 



(2.23) 
(2.24) 

(2.25) 
(2.26) 
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wher e Erfc is the complementary error function. To get the parallel integral in In 
[Eq. (2.14)], we simply take the limit of / as a — > 1. We then obtain the following: 



POO 

In- = / dx± x± exp[— x] 
Jo 



7T 



exp[— oj\ 



Sgn[Im[uJN)]i + Erfi[vu; 



(2.27) 
(2.28) 

Each of the terms in Eqs. (2.27) and (2.28) can be evaluated in a straightfor- 
ward manner (using a handbook of integrals or a symbolic integration package, for 
instance). The resulting equations are 



f°° 71 

In,+ — / dx± x± exp[— x\] expjo;] — ^=Erfc[Vo)]. 



l N, 



— V / 7r^2exp[— Cj] ( Erfc[V— 



7r 3 2 exp[u>] ( Erfc [v Cj 



(2.29) 
(2.30) 



where Cj = —ujn/C,- Plugging these expressions into the original integrals of interest 
from Eq. (2.12), we get 

uon + 1 



cf v r-Fjtfe = n ae - (uj n + 1) 

u} N + \L n /R\(xl + a?j2) 2 



R 



L, 



exp[— Cj] ^Erfc[ 



-co 



(2.31; 



for electrons and 

7-V U N ~T 



uj N - \LJR\t (x\ +x 2 ± /2j 



Mi 



n 0i n 
' r 2 



(uj n - r) 



R 



exp[Cj] I Erfc [v Cj 



(2.32) 

for ions, with r = To«/To e the ratio of ion to electron temperatures. Substitut- 
ing Eqs. (2.31) and (2.32) into the quasineutrality expression (2.12) results in the 

R 



following dispersion relation: 



1 + r- (r - uj n ) 
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- r (1 + u) N ) 



71 



exp[u>] ( Erfc [v a) 
R 

exp[— Cj] (Erfc[ 



-0J 



(2.33) 



0. 



The presence of the complementary error functions in the above dispersion 
relation makes analysis complicated. We simplify matters by assuming \Cj\ <C 1 (i.e. 
\ujn\ <C |£|) and taking r = 1. To lowest order in u N , we find 



n\R/L n \-2 , „ 

'UN = 7= 3/2 I 1 - V-lj 



'tor\R/L n \ 



(2.34) 



This expression must be treated carefully. Depending on the choice of branch cut, 
\J — 1 = ±z. The choice of branch cut, coupled with an assumption about the sign of 
Im[u N ], also sets a restriction on the signs of the real and imaginary parts of Jujn- 
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First, we take a branch cut along the negative real axis, so that the arguments 
of complex numbers are defined on the interval [— it, it). In this case, y/—l = —i and 



'UJ N 



7i\R/L n \-2 
\R/L n 



3/2 



(2.35) 



If Im[ujff] > 0, then Re^LO^] > and Im[y/Ufj\ > 0. Consequently, Eq. (2.35) is 
only valid when \R/L n \ > 2/n. If we instead assume that Im\u^\ < 0, then we have 
^uj~n] > and Imly/uJ^] < 0. This is clearly not possible in Eq. (2.35) since 
/aJjv] = lTn[- s /uJ^\, so no damped waves exist for this choice of branch cut. 



Re[ 
Re\ 



Now we take our branch cut along the positive real axis, so that the arguments 
of complex numbers are defined on the interval [0, 2ir). For this case, = i, and 
the equation for tJuJn becomes 



'UN 



n \R/L r 



8ir\R/L n 



3/2 



(2.36) 



If JmfcuTv] > 0, then Re[y/ujJJ} > and Jm[ v /cJ/v] > 0. Since Re^uj^] = — Irn[^/uJ^\ 



in Eq. (2.36), no growing modes are allowed for this choice of branch cut. For 
ira[a>jv] < 0, we get Re[y/uJ N ] < and Im[yJu)N\ > 0. This is only satisfied for 
\R/L n \ < 2/tt. 



Combining the results of Eqs. (2.35) and (2.36) and keeping in mind their 



range of validity, we obtain our solution for u^'- 



(tt \R/L n 


-2) 


TT \R/L n 


- 2 | 


4vr \R/L n 


3 



(2.37) 



which indicates instability for gradients steeper than the critical gradient, given 
by \R/L n \ crit = 2/tt. A similar, if somewhat messier, calculation can be done for 
arbitrary temperature ratio and with lowest order finite Larmor radius effects in- 
cluded |38] . The result is 



oj- 



[(1 + r) (| \R/L n \ - 1) - kip! \R/L n \ (tt/2 " 1)] ! 



VthA 



2tt (1 + r 3 ) 2 \R/L t 



(r 2 - 1 ± 2r 3 / 2 ^) k ± (H, 
(2.38) 

where the + sign applies for \R/L n \ > \R/L n \ crit and the — sign applies for \R/L n \ < 

\ R / L n\ cri V Wllere 

\ R I L ^U = 7r( l + r) _ A;i p 2(7r _2)- (2 ' 39) 



2.3 Comparison with fluid theory 

Previous studies [121 1131 EH ESI ESI EI] of the low-/? Z-pinch system considered 
here have identified two distinct linearly unstable modes. The first mode is the well- 
known ideal interchange mode described by MHD, which is the dominant instability 
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(with growth rate scaling like 7 ~ c s /\/RL n , c s = a/ (T 0e + T 0i )/mi) in the strong 
gradient regime. However, the ideal interchange mode becomes stable at moderate 
gradients, leaving only the shorter wavelength entropy mode. This mode consists 
of perturbations to both the plasma density and temperature, but not the plasma 
pressure. 

In Fig. 2^ the dependence of the growth rate on the density gradient scale 
length is considered for a number of different models. Collisional and collisionless 
MHD equations accurately predict the growth rate in the regime of strong density 
gradient, but they do not describe the dynamics of the entropy mode and thus 
are not applicable for sufficiently weak gradients (where the interchange mode is 
stabilized). Fluid models capture the correct qualitative behavior for the entropy 
mode, but they do not provide good quantitative agreement with the growth rate 
and they underestimate the critical gradient determined from the gyrokinetic model 
by a factor of approximately two. 

We see in Fig. 2J2 why a fluid description is not necessarily sufficient for this 
system: the distribution of particles in velocity space possesses nontrivial structure. 
In particular, there is an energy resonance in the velocity space arising from the 
presence of curvature drifts. Additionally, fine scale structure is present in the 
distribution of perpendicular velocities. The development of such fine scale structure 
is discussed in more detail in Chapter |4j 

As a final note, we point out that even for this highly simplified system, the 
kinetic stability calculation was quite involved. To accurately model linear physics 
for kinetic systems of practical interest to the fusion community, we would need to 
retain finite Larmor radius effects, finite temperature gradients, trapped particles, 
magnetic shear, and a host of other important features. Such calculations are bur- 
densome at best and quite often analytically intractable. Furthermore, while linear 
analysis is undoubtedly useful in providing qualitative insight into the gross charac- 
teristics of plasma behavior, a nonlinear treatment is of course necessary to address 
turbulent dynamics. The inherent difficulties of analytic calculations make com- 
puter simulations a critical component in advancing our understanding of plasma 
turbulence. 
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Figure 2.1: Growth rate of the interchange and the entropy mode as a function of 
L n /R for two different values of kp s (both with r = 1). The collisionless gyrokinetic 
growth rate (red " x " marks) is compared to the growth rates from the gyrofluid 
model (black dotted-dashed line), the ideal collisional (green dashed line) and col- 
lisionless (gree solid line) interchange mode, and the collisional (blue dashed line) 
and collisionless (blue solid line) fluid entropy mode. The kinetic model is necessary 
to obtain the correct stability boundary and to obtain the correct growth rate for 
weak to moderately strong gradients. Figures taken from Ref. [55] . 
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Figure 2.2: Real (a,b) and imaginary (c,d) parts of the ion (a,c) and electron (b,d) 
velocity distribution functions for the case of a moderate density gradient (L n /R = 
0.5) and large kp s (= 38). Axes are normalized to Vth, s - We see significant structure 
in v± for ions in addition to an energy resonance arising from the curvature drift. 
Figures taken from Ref. [38] . 
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Chapter 3 

Turbulent transport hierarchy for gyrokinetics 



3.1 Introduction 

As noted in Chapter [TJ a full treatment of turbulent transport in weakly col- 
lisional plasmas involves the interaction of slowly evolving, large scale profiles and 
rapidly evolving, small scale turbulence. It is not straightforward to neglect either 
process, because the large scale profiles and small scale turbulence are dynamically 
coupled. The evolution of the profiles is governed by the turbulence, and the in- 
stabilities that give rise to the turbulence respond sensitively to the profiles. The 
parameter space that characterizes these interactions is very large, limiting the ulti- 
mate applicability of parametric fits (271 HI] . The analytic and numerical difficulties 
associated with this wide range of time and space scales are compounded by the 
kinetic nature of the instabilities driving the turbulence; in principle, one must con- 
sider a six-dimensional phase space consisting of three dimensions each for physical 
and velocity space. Such a system is both analytically and numerically intractable. 
Consequently, it is necessary to work with a reduced model. 

In this chapter, we will derive such a reduced model by taking advantage of the 
wide space and time scale separations present in many weakly-collisional plasmas. 
In particular, we closely follow the treatment of Refs. (Ml 06] to introduce a set of 
ordering assumptions (the 5f gyrokinetic ordering [SI [9], [10]) that leads to a reduction 
of the phase space and to the development of a hierarchical set of equations, in which 
equations describing the turbulent fluctuations and equilibrium profiles are coupled, 
but evolved separately. This represents a significant simplification of the system, 
which we show in later chapters allows for computationally efficient, first-principles 
simulations of turbulent transport evolution over long time scales. 

3.2 Ordering assumptions 

We take the Fokker-Planck equation as the starting point from whence we will 
derive our hierarchical set of equations describing turbulent transport: 

where f s represents the distribution of particles of species s in position r and velocity 
v, q s is particle charge, m s is particle mass, E and B are the electric and magnetic 
fields, c is the speed of light, and C[f s , f u ] is a bilinear operator describing the effect 
on particles of species s of collisions with particles of species u. For convenience 
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of notation, we shall henceforth drop the species subscript s where it leads to no 
ambiguity. 

Assuming the potential energy of nearest-neighbor interaction is much less 
than the kinetic energy of the particles, this equation describes the full range of 
dynamics in a six-dimensional phase space and time for particles of species s moving 
in a self-consistent electromagnetic field. It does not take into account the effect 
of external sources of particles, momentum, energy, etc. These will certainly be 
present in fusion devices, but they are local processes, and their time and space 
scales are associated with thermodynamic processes. Consequently, we neglect them 
in our analysis and insert them ad hoc when we consider the slow evolution of 
thermodynamic quantities such as density and temperature. 

Unfortunately, the comprehensive nature of the Fokker-Planck equation makes 
theoretical analysis burdensome and numerical simulations computationally infeasi- 
ble. To make progress, we simplify our model by adopting a variant (HI US] of the 
5f gyrokinetic ordering [8j El HO], which exploits scale separation in time and space. 
We first separate quantities into equilibrium and fluctuating parts: 

f = F + 5f, B = B + 5B, E = 5E, (3.2) 

with 

-L „ „ e « i. (3.3) 

-TO -DO 

Formally, we define the smallness parameter e as a ratio of small to large spatial 
scales within the plasma: 

e = £, (3.4) 

where L is a typical scale length associated with the equilibrium and p is the radius 
of particle gyration in the equilibrium magnetic field. We separate spatial scales 
by assuming cross-field fluctuations vary on the gyroradius scale, while all other 
quantities vary on the equilibrium scale: 

V ± 6f ~ ^ (3.5) 
P 

VF ~^, V||<5/~^. (3.6) 

We separate temporal scales by assuming gyromotion is faster than the dynamic 
frequencies of interest, which are themselves much faster than the evolution of the 
equilibrium profiles: 

7-- 1 ^ e 2 uj ~ e 3 n 0i , (3.7) 
where f2 i = \q\B /rriiC is the ion gyrofrequency and 

dF F 38 f 



uSf. (3. 



dt t ' dt 

Additionally, we take the the Ex B velocity to be an order smaller than the thermal 
velocity: 

c5E 

evth, (3.9) 



B 
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where we have assumed no background electric field (and therefore no equilibrium 
flow). We can express the fields in terms of potentials as follows: 

1 <9 A 

E = -V$--— (3.10) 
c at 

B = VxA (3.11) 

5B = Vx<5A (3.12) 

where $ is the electrostatic potential and A is the vector potential. Note that our 
ordering requires 5A/A ~ e 2 (since VA ~ A /L, V5A ~ SA/p, and 5B/B ~ e), 
so that the electric field is primarily electrostatic in nature. 

The final ordering assumptions we make are that the collision frequency v is 
comparable to the fluctuation frequency and that the distribution function varies in 
velocity space on the scale of the thermal velocity: 

9fs fs / O 1 O \ 

V ~ U, -7— ~ • (3.13) 

<9v v th)S 

Note that this choice of ordering does not prevent us from considering the cases 
<Cw and i/>w as subsidiary orderings [H] . In general, subsidiary orderings can 
be applied using a large number of plasma parameters as the expansion parameter. 
While most of the potential expansion parameters have order unity variations across 
many experiments, some, such as the electron-ion mass ratio, the plasma beta, and 
the aspect ratio of the device, are good expansion parameters for a large class of 
systems. However, we are interested in deriving a set of equations that are widely 
applicable to turbulent transport studies, so we do not consider such subsidiary 
expansions, which would limit the range of validity of our model. 

With these ordering assumptions, we can now expand the Fokker-Planck equa- 



tion (3.1) in the smallness parameter e. What we find is a set of ordered equations 
that ultimately provides us with information about the evolution of the equilibrium, 
instabilities, fluctuations, and transport. In what follows, the ordering is taken rel- 
ative to v t hF /L. 

3.2.1 Fq does not depend on gyrophase 

The lowest-order equation in our e expansion is 

-vxB o .f°=-^ = 0, (3.14) 
mc ov ov 

where $ denotes the gyroangle. Thus, the equilibrium distribution function is in- 
dependent of gyroangle. We note that in obtaining the first equality in the above 
equation, it is convenient to use cylindrical (v±, v\\) or spherical £) coor- 
dinates in velocity space (£ = v\\/v is the pitch-angle). Both of these coordinate 
systems will be used frequently in this and later chapters. 
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3.2.2 Fq is Maxwellian and 5f can be decomposed use- 
fully 

At the next order (e°), we find 

q f v x BA 9F vx B d/i „ , 1K v 

v-VF +v r V/i + - Ei H •— — H — = > C Fo,4 , (3.15) 

m V c J ov c ov ^— ' 

where we are using the notation Sf = /1+/2 + ..., etc., with f n /F ~ e n . Multiplying 
this equation by 1 + In F and manipulating, we have 



V • (F lnF ) v + (1 + lnF ) v ± ■ V/i + ^ 
= (1 + lnFo 



* ( El + ^]F lnF 

m \ c 



(3.16) 



(^ + £C[F ,F .jj. 
Integrating this expression over all velocities gives 

J d 3 v ( V • (F In F ) v + (1 + In F ) v ± • V/i = £ In F C[F , F 0u ] J , (3.17) 



where we have used the divergence theorem to eliminate the third term in Eq. (3.16 ), 
and we have used the fact that collisions locally conserve particle number to assert 
Jd 3 vC[F ,F 0u } = 0. 

Before proceeding, we need to define an intermediate spatial average, which we 
will see is equivalent to a flux-surface average. We restrict our attention to axisym- 



metric equilibrium magnetic field configurations (Fig. 3.1 ), which can be represented 
in the following general form (a fuller discussion of the magnetic geometry is given 
in Appendix [A]): 

B = I (if)) V0 + W> x V0, (3.18) 

where if) = (27r)~ 2 J dVB ■ V# is the poloidal flux, dV is the volume element, and 
9 and are the physical poloidal and toroidal angles respectively. The quantity 
I(if>) — RBt is a measure of the toroidal magnetic field, where R is the major 
radius and B T is the toroidal magnetic field strength. 

The intermediate spatial average of any quantity J-(y) is then defined as fol- 
lows: 

1 p2tt pit pipQ+A-ip/2 

«.F(r))) = - / # / d6 / # ^(r), (3.19) 

V JO J-n Jip -Af/2 

where 

p2n pit fipo+Aip/2 

V = d<f> d9 dip J (3.20) 

JO J —n J^ Q -A.ip/2 

J = (Vtjj x V0 ■ V0)" 1 , (3.21) 

and ipo denotes location in the radial coordinate if). Formally, we define the in- 
termediate spatial length Aip to be order e 1 ' 2 L (i.e. p <C Aif) <C L). Note that 
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Figure 3.1: Schematic of an axisymmetric magnetic field configuration. The flux 
surface, labeled by pressure p or toroidal/poloidal flux \l/ has no variation in the 
toroidal (ip) direction. Figure taken from Ref. |47j . 

in the limit where the volume V becomes vanishingly small, the spatial average of 



Eq. (3.19) reduces to the usual flux surface average. We also point out that the 
spatial averaging is done about a fixed point in space so that it does not depend on 
time. It is possible to define the spatial average so that it is taken with respect to a 
fixed flux surface. However, since the flux surfaces themselve evolve in time on the 
equilibrium time scale, this would require additional terms to account for the time 
dependence. 



We now apply the intermediate spatial average to Eq. (3.17). The divergence 
theorem applied to the first term gives 

r d\ 

((V-(FolnF Q )v)) = J — -vFolnFo, (3.22) 

where dA is the area element whose vector direction is normal to the surface bound- 
ing the volume integral (the flux surface). Since the magnetic field lies within the 
flux surface, dA • Bo = 0, so that 

d\ / — • vF In F = / d 3 v / — • v ± F In F = 0. (3.23) 

The last equality follows from the fact that Fq is independent of gyroangle, and vi 
is an odd function of gyroangle. 



The second term in Eq. (3.17) can be integrated by parts: the surface term 



is zero at this order and the remaining term is dropped in ordering by e 1//2 . Conse- 
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quently, we find 



J d 3 v InFo ^C[F ,F( 




0. 



(3.24) 



From Boltzmann's if -Theorem, the only solution to this equation is F = Fm, where 
Fm is a Maxwellian in velocity. 



Plugging F = Fm into Eq. (3.15), we find 

dfi 



v± ■ v/i - fi 



v • VF, 



M 



v ± • V 



where T = mvf h /2 is the equilibrium temperature, and the v\\V- ^ — 



t 



(3.25) 



Fm term has 



been neglected at this order. The velocity space derivatives thus far have implicitly 
been taken at fixed particle position r. However, it is now useful to switch to the 
guiding center variable R = r — p, where p — b X v/Oq is the gyroradius vector: 



d_ 



d_ 



R 



dR 



d_ 
dR 



d_ 



R 



(3.26) 



where subscripts on the derivatives denote quantities that are held constant during 
the differentiation. It should be noted that there is no ambiguity in the use of the 
V notation, since d/dr = <9/<9R. Using this result in Eq. (3.25), we have 



-fir 



dfi 



-v ■ VF M - Yi ■ V 



R 



To 



(3.27) 



The homogeneous solution h satisfies 



dh\ 



(3.28) 



telling us that h is independent of gyroangle at fixed guiding center position. 

We next proceed to find the particular solution. Applying Eq. (3.26) to the 
righthand side of Eq. (3.27) and remembering that F Q is independent of gyroangle, 
we find 



dfx 



fin 



R 



d_ 



R 



M 



Upon gyro-averaging, we obtain 



b - VF, 



M 



1 + 



0. 



v\\h-VF A/l 



M ■ 



(3.29) 



(3.30) 



which is simply a statement that the equilibrium distribution function is constant 



on a flux surface. Using this result in Eq. (3.25), we have 

dfi 



v± ■ V/i - fi 



-v± ■ V 



M 



1 + 



(3.31) 
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which has the particular solution 

fi P = ~F M -p-VF M . (3.32) 
Jo 

If we redefine our equilibrium distribution function by absorbing this term, we see 
that we have the beginning of a Maxwell-Boltzmann distribution for guiding centers: 

F = F M (R)exp^-^ . (3.33) 

Our solution then has the form 

f = F + h + f 2 + ... (3.34) 
3.2.3 The gyrokinetic equation 

At order e 1 , we derive the gyrokinetic equation. It is convenient to transform 
to gyrokinetic variables, which is easier if we start again with the Fokker-Planck 
equation: 

Of dR df dedf dfidf dfidf _sr^ n \ f fl /o ocr\ 

m + ~dt ' M + Tt~d~e + dt + dtM - 2>l/' ™> V-W 

u 

where e = mv 2 /2 + g$ is the particle energy and // = mv\/2B is the magnetic 
moment. The (9(e) terms yield 

dh dR d F de df 2 r , , , 

m + ^-M {Fo + h) 'T di- n °M =c[h + p - VF ^ (3 - 36) 

where we define 

C[f] = (C[/, + C[F M , /„]) . (3.37) 
We take the gyroaverage of this equation (at constant R) to eliminate the f 2 term: 
| + (-) r . V(Fo + A) _|(|)^ (C[A])r , ( ,38) 

where we have used the fact that both F M and h are independent of gyroangle at 
fixed guiding center position. To the order we will need, the guiding center velocity 
and power are 

= ^l|b+(v D ) R (3.39) 

\ uv I R 

and 



\ dt i R 



where 



b 
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(3.41) 



is the gyroaveraged guiding center drift velocity. Substituting Eqs. (3.39) and (3.40) 

qF Q fd(x) R 



in Eq. (3.38), we obtain 
dh 



, )f +V]l b.Vh + (v D ) R .V(F + h) = (C[h]) K + 



dt 



v 9Ac 
c ' dt 



. (3.42) 



To obtain the gyrokinetic equation in a standard form, we first split h into two 
pieces: one that varies on the equilibrium spatial scale (called the neoclassical part) 
and one that varies on the gyroradius scale (called the turbulent part): 



h = ht + h r 



(3.43) 



We then apply the intermediate spatial average defined previously to obtain an 
equation for the evolution of h nc : 



dhnc , u Y7h , up in\h 1\ qF ° V dA ° 
+ v\\h- V/inc + v B • VF = {C[h nc \} R - 



dt 



Tn c dt 



where 



vlh ■ Vb + 



2 B 



(3.44) 



(3.45) 



consists of the sum of the curvature and V-B drifts. Subtracting the neoclassical 
equation (3.44) from Eq. (3.42) yields the gyrokinetic equation for the evolution of 
the turbulent distribution function: 

^ + vfi ■ Vh + (v x ) R • V(F + h t ) + v B ■ Vh t = (C[h t ]) n + (3-46) 



with 

c 

X B 

being the generalized E x B drift velocity. 



b x Vx 



(3.47) 



3.2.4 Transport equations; thermodynamics 

The 0(e 2 ) equation includes terms involving f'2, for which we have no expres- 
sion. The goal in this subsection is to manipulate the equation to eliminate second 
order quantities in favor of products of first order quantities. We will accomplish 
this by taking a moment approach and averaging over intermediate space and time 
scales. The result, as we will see, is a closed set of fluid equations for the evolution 
of the equilibrium density and temperature profiles. Because of the time and space 
averaging, these equations are only valid when there are no important time or space 
scales between the turbulence and equilibrium scales. 

Before we begin the moment approach, we define the intermediate time average 
of a quantity T(t) in a manner analogous to the intermediate spatial average of 



Eq. (3.19): 



I r T °+— 



JF( T0 ) = — / dtF(t), (3.48) 



where u 1 <C At <C t, and we formally define At ~ er. 
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3.2.4.1 Slow density profile evolution 

To obtain an equation for the evolution of the equilibrium density profiles, we 
first take the density moment of the Fokker-Planck equation. The term involving 
the Lorentz force vanishes because it is a perfect divergence in velocity space, and 
the collisional term vanishes because of local particle conservation. We are left with 
the usual continuity equation: 



d 3 v 



df 
dt 



V • (v/) 



0. 



(3.49) 



Applying the gyrokinetic ordering from the beginning of the chapter, we see that 
both 8/2/ dt and V • (v/ 2 ) enter at order e 2 . Since we do not want to solve for / 2 , 
we must eliminate it from the equation. We will accomplish this by averaging over 
intermediate scales in space and time. 

Performing the intermediate spatial average of Eq. (3.19) and using the diver- 
gence theorem, we can write 



«V-(v/)» 



dA 



v/ 



dA Vip 



v/. 



(3.50) 



By inspection, one can see that this can in turn be written 



dA W 



1 d 



(3.51) 



V |W>| J Vdi\) 

We can manipulate this into a more useful form in a couple of steps. First, we use 
the expression for the axisymmetric magnetic field given in Eq. (3.18) to get 

(v • VV) / = -R 2 V<P ■ (v x B ) /. (3.52) 

Noting that v x Bq = — ^o§| and integrating the righthand side of the above 



expression by parts in gyroangle, we obtain 

J d\ (v • W) / = J d\ (R 2 V(p ■ v) 



(v x B 







dj 
<9v 



(3.53) 



Using Eqs. (|3.50|), (|3.5lD, and (|3.53|) in Eq. ^3.49|), we obtain 
d\ . ) } 



dt 



1 d 



Vdip 



V J d\ (R 2 V<p ■ v) 



v x B 



0J 



dl 
9v 



0. (3.54) 



Substituting for (v x Bo) ■ df/dv from the Fokker-Planck equation (3.1), we have 

Vd~i) 
Q 



,3 df 



+- 



m 



E + 



v x SB 



dl 
dv 



(3.55) 



0. 
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We now apply the intermediate time average defined by Eq. (3.48), plug in our 
expression (3.34) for /, and examine each of the resulting terms up to order e 2 . 



The first term in Eq. (3.55) becomes 



df 

dt 



dFo dh dh 
dt dt dt 



(3.56) 



Since h and f% both vary on the fluctuation time scale, the only term that survives 
the time average is the one involving F : 



d 3 v 



df 
dt 



dn 
dt ' 



(3.57) 



where = ((f d 3 vFo)). Because of the slow time variation of Fq, this term enters 
at 0(e 2 ). The other term in Eq. (3.55) involving df /dt is treated analogously. 



However, the prefactor multiplying it is of order e, dropping the overall order to e 3 . 
Consequently, we may neglect it. 



We next treat the term containing v ■ V/. Employing Eq. (3.34) for /, we 



have 



Vd~i) 

" vd~i) 



v 



d 3 v— (i? 2 V</>-v) (v-V/) 



V 



777 C 

d\— (i? 2 V0 • v) (v • V [Fq + h + f 2 }) 



(3.58) 



Since Fq = Fm(iP) and • V"0 = 0, the term in the integrand containing Fq is 
odd in velocity space and therefore integrates to zero. Integrating the f'2 term by 
parts in space, we see that the gradient operator is transferred to the equilibrium 
quantity R. Due to the slow cross-field spatial variation of equilibrium quantities, 
this drops the term to C(e 3 ), so it does not contribute at this order. 

We can also integrate the h term by parts in space, and we find that the 
turbulent piece vanishes due to periodicity. All that is left is a term involving the 
neoclassical piece of h, which we address now. First we define the neoclassical 
pressure tensor: 



d wvL 



d s v 



w 2 bb + 



(I - bb) 



h + P 



(3.59) 



where I is the identity tensor and P = f d 3 v (vyvj. + Vj_V||) h nc is an antisym- 
metric tensor containing the off-diagonal components of the pressure tensor. The 



neoclassical term from Eq. (3.58) can then be written 



1 d 



vm\ V I d3v f^- v)(v - Vh 

1 mc d I 
V q dip\ 



' ""' ' ((V [V • (P • R 2 V<£) - P : V (i? 2 V0)] 



(3.60) 
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where the double dot tensor product is defined 



(3.61) 



The first term on the right in Eq. (3.60) vanishes upon application of the divergence 



theorem due to periodicity in toroidal and poloidal angles. Because V (RV(j)) = 
(RV0 — V0R) is an antisymmetric tensor, the symmetric part of the pressure tensor 
vanishes in the second term on the right, leaving 



IA 

Vdi) 



((y J d\™ (i? 2 v0-v) (v-vh m 

P : V (i? 2 V0) 



1 mc d l . ^ . 
V q dip\ 



(3.62) 



0. 



To obtain the final equality, we made use of the identity J d 3 R J R d 3 v = J d 3 r J r <i 3 v, 
where the subscripts on the integrals denote the variable to be held fixed. Since h nc 
is independent of gyroangle at fixed R and (v|jVj_ + Vj_vii) is odd in gyroangle, we 
see that the P inside the spatial average of Eq. (3.62) vanishes. 



Now we examine the term containing the Lorentz force: 



IA 



{(v j d\ (R 2 V(f) ■ v) (c5E + v x SB) ■ ^ 



IA/ 



V [ d\ (R 2 \7(f) ■ v) (cSE + vx5B)--^-(F + h + 5f 2 ) 
J ov 



(3.63) 



First, we consider the terms containing F . For the magnetic force term we find 

J c/ 3 v (R 2 V(f) -v) (vxffi) — ~ J d\ (R 2 V(j) • v) (v x SB) ■ vF = 0. (3.64) 

For the electric force term, we have 

dF 



d 3 v (R 2 V(f) ■ v) SE ■ ^ = f d\ (R 2 V(j) • v) A . f F SE \ 
o\ J ' av 



d\R 2 V(j) 



V$ + -^(A + (5A) 

c ot 



(3.65) 



where we have integrated by parts and substituted Eq. (3.10) for the electric field 



between lines one and two. Upon application of the intermediate time average, the 
SA term drops in ordering. We next apply the intermediate spatial average: 



R 2 V<p ■ V$F 



^0 



0. 



(3.66) 



where the final equality follows from the fact that the equilibrium is axisymmetric 
and all other quantities are periodic in torodial angle 0. 
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Now we consider the term involving h. Integrating by parts in velocity space, 



we get 



d 3 v (R 2 V(f) ■ v) ( 5E 



v x 5B 



c 



*-/*(*V,..)». 

d\R 2 V<p- [SE + 



5E + 

c 



vx5B 



h. 

(3.67) 



h 



Rewriting 5E and <5B in terms of the potentials $ and A, we have 



5E + 



v x SB 



^ IdA v 

-V$ — + - x (V x <JA) 

C OT C 

i BA 1 

-V$ - - V + " [V (v • 6 A) - (v • V) <5A] 

c at c 

„ 1 <9A /v \ . . 



(3.68) 



where x is defined in Eq. (4.5). The terms involving v»h ■ V5A and dA/dt are 



higher order and can thus be neglected (the dA/dt term was retained earlier when 
multiplied by F , but here it is multiplied by h). Further, only the fluctuating part 



of h contributes, since the intermediate spatial average of Eq. (3.67) drops the order 



of the h nc part. The • V5A term vanishes by integrating by parts in space and 



using the identity from Eq. (3.26): 



J d 3 v R 2 V(f) ■ [(vl • V) 6 A) = J d 3 v R 2 V(p ■ 5 A (v ± ■ Vh) 



0. 



where we have used the fact that h is independent of gyroangle at fixed R and 5 A 
is independent of gyroangle at fixed r. 



is 



ependent of gyroangle at fixed r. 

The only part of the distribution function contributing to the collisional term 



h"nc — h"nc P ' F'l 



(3.70) 

because h t vanishes upon spatial average and C[F M ] = 0. Collecting results, we 
have the equation for the evolution of the equilibrium density: 

3.2.4.2 Slow temperature profile evolution 



(3.71) 



The equation for the evolution of the equilibrium temperature profiles is de- 
rived in a manner analagous to the equation for the equilibrium density profiles. 
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We begin by taking the energy moment of the Fokker-Planck equation (3.1) and 
performing the intermediate spatial average of Eq. (3.19): 



cfv 



mv 



!+v.V/+« 
at m 



E 



v x B 



d_l 
<9v 



(3.72) 



We then substitute Eq. (3.34) for / into Eq. (3.72), average over the intermediate 
time scale, and consider the order of each of the terms. 

As before with density evolution, the terms involving time derivatives of h and 
5/2 do not contribute at this order. The contribution to the first term of Eq. (3.72) 
from F is 

3 dn T 



,0 mv 2 df 

cfv 

2 dt 



2 dt 



(3.73) 



where T = ((/ cfv mv 2 Fo)) /2n . 

The term in Eq. (3.72) containing v • V/ can be manipulated as it was in 
Eqs. ( |3~50| -( |3~53| to obtain 



mv 



cfv —v ■ V/ 



V84) 



V d\ - ° 2 



R z V<j> ■ v 



x B ) 



df 

<9v 



(3.74) 

Substituting for (v x Bo) • df /d\ from the Fokker-Planck equation (3.1), we have 



cfv 



mv 



v V/ 



V / c/ 3 v 



mv mc 



+- 



m 



E + 



vx5B 



c 



9/ 
9v 



(3.75) 



Using the same methods employed to derive the density evolution equation, one can 
show that the terms involving df /dt and v • V/ do not contribute at this order. 

We focus first on the Lorentz force term. As before with density evolution, we 
use / = F + h + 5 $2 and examine each term. First, we consider the terms involving 
Fq. The magnetic force term is zero since dF /dv ~ vF and (v x <5B) • v = 0. The 
electric force term gives 



J cfv (i? 2 V0 • v) ^c5E 



OFn 



<9v 



rf 3 vci? 2 V0 



mv 



5E + (v • SE) v 



F 



(3.76) 



where we have integrated by parts in velocity space. Writing the electric field in 
terms of potentials and noting that the intermediate spatial average of fluctuating 
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quantities is zero, we find 



J rf 3 v (i? 2 V0 ■ v) ^c5E • 



dFo 
<9v 



dt 



(3.77) 



Next we consider the term involving h: 
J d\ (# 2 V0 • v) ^ (c5E + v x 5B) 



9v/ 



d\R 2 V<p 



mv 



(c5E + v x 5B) + mcv (v • 5E) 



(3.78) 



where we have integrated by parts in velocity space and neglected the h nc terms 
because the spatial average makes them higher order. Also, only the electrostatic 
part of the electric field contributes at this order. Rewriting the above expression 
in terms of potentials, we have 

,2 



cfvirv</> 



mv 



(cVx + v_l • V±5A) — — h mcv (v • V$) 



ht 



(3.79) 



where the term involving v»b ■ V5A has been dropped since it is higher order. The 
• V5A term can be neglected following the same argument given after Eq. (3.68). 
The part of the collision operator that contributes at 0(e 2 ) in Eq. (3.75) is the same 
as for particle transport, C[h nc \. 

Next, we consider the Lorentz force term in Eq. (3.72). The magnetic force 
term can be written as a perfect divergence in velocity space and subsequently 
vanishes upon integration. However, due to the presence of the v 2 factor, the electric 
force term is nonzero. Integrating this term by parts in velocity gives 



mv 2 q „ df 
2 m av 

Considering the electrostatic part of E, we have 



d 3 v qE ■ v/. 



(3.80) 



q / d 3 v (v ■ V$) f = q d\ 



q / d 3 v 



d(&f) 9$ 

v.v(*/) + -L£- % -* 
1 J ' dt J dt 



df 

dt 



+ v-V/ 



(3.81) 



where we used the continuity equation (3.49) to obtain the final equality. The 
intermediate time average eliminates the second term and the Fo piece of the final 
term on the last line. Performing the intermediate spatial average eliminates the 
neoclassical terms and the F part of the first term, leaving 



q 



d 3 v (v • V$) / 



q 



v • V ($ht 



9$ 

~dt 



(3.82) 
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We can place the first term in a more convenient form by first using Eqs (3.50 )-(3.53 ): 

1 d 



J d 3 v q\ ■ V {$ht 



V J d 3 v (i? 2 V0 • v) g$ 



v x B 



Noting that (v x B ) ■ dh t /dv = —B (dh t /d'&) r and using Eq. ( 3.26[ ), we get 



(3.83) 



J d 3 v gv • V ($h t ] 



Sin 



(i? 2 V0 • v) v • Vh t 



(3.84) 



where we have used (dh t /d$)n = 0. Integrating by parts in space then gives 
d\ ^^(^ t ))) = -~ 



V J d\ {R 2 V(p ■ v) mc h t v ■ 



(3.85) 

which cancels with the last term in Eq. (3.79). 

We now consider the inductive part of the electric field, Ej. Since F M is 
isotropic in velocity space, we have 

q& 



d A w q (v • Ej) F M ( 1 — 



0. 



The term involving f 2 does not enter at this order, so we are left with 



d 3 v q (v ■ Ej) (p • VF M -h)= / rf 3 v - I v 



OA 

c V 

Applying the intermediate spatial average, we obtain 



(h-p- VF 2 



Mi 



(3.86) 



(3.87) 



J i^jIvEj) (p-VF M -h) 



3,. 9 



MA 



^+ V 



5Ao 

at 



(3.88) 



Combining this expression with the one remaining term from the electrostatic part 



(i.e. the last term in Eq. (3.82)) gives 



3.- 9 



~m ht 



at 



h, 



(3.89) 



The final term to evaluate in Eq. (3.72) is the inter-species collisional energy 
e following standard form (see e.g. Ref. [67]): 

^v^£C[/,/ u ]\ \ =J2^: u (T 0u -T ), (3.90) 



where 



{m s m u ) l/2 q 2 s q 2 u n u In A 



v° u = 6.88 



(m s T u + m u T s 



,3/2 



1.55z/ s 



(3.91) 
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with 



4nn u q^q; In A 



(3.92) 



rnj 2 (2T s f /2 

the collision frequency, q the particle charge, and lnA stJ the Coulomb logarithm. 
Collecting results, we have an equation for the evolution of the equilibrium pressure, 
p Q = n T : 



3 dp 1 d 
2~dt = Vdi) 



V / d\ 



mv 



d v q 



-R 2 V(p 



' dx dA ~ 

ht di- v --m hn 



^Fm + h t Vx - —C[h nc ] 
at q 



+ 5>o^r (T 0u -T 



(3.93) 



3.2.4.3 Species-summed pressure equation and turbulent heating 



The term ((J <i 3 v q s h t>s {dxi 'dt))) m Eq. (3.93) describes turbulent heating of 



the equilibrium. In order to illuminate the nature of this turbulent heating, we sum 



Eq. (3.93) over species and consider the evolution of the total pressure, px = ^2 s Po s '- 



1 d 



y i)t • 

s 



3 m s v 2 R 2 m s 
V I a v 



q s 



{V<f>-v)C[h r , 



dx <9A ~ 
at at 



(3.94) 



where 



(3.95) 



is the species-summed turbulent heat flux. We now proceed to show that species- 
summed turbulent heating term, ((f d 3 v q s h tjS (dx/dt))) is zero for steady-state 
turbulence. 

First, we consider the electrostatic component of X'- 



<9$ 

~dt 



q s J d\ ht 



q 2 <9$ 2 
2T 0s dt 



(3.96) 



where we have used quasineutrality. Upon time averaging, this term vanishes. Next 
we consider the inductive component of x : 



105 A 



c dt 



£ q s J d 3 v vh ty 



185 A 



■53 



c dt 
-((5Ej-53)) 



(3.97) 
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where 53 is the perturbed current, and 5~Ei is the inductive part of the fluctuating 



electric field. Applying Ampere's Law and using the Coulomb gauge, Eq. (3.97) 
becomes 



ld5 A 



c dt 



GTV v/l 



t,s 



1 85 A 



c dt 

ld_ 
2di 



V 2 5A 



(V5A : (V5A) T ) 



(3.98) 



where the superscript T denotes the matrix transpose. This term also vanishes upon 
time averaging. Consequently, we have 



E 



d vg s /i tjS — 
dt 



0. 



(3.99) 



where the usual intermediate time average is implied. Therefore there is no species- 



summed turbulent heating. Using Eq. (|3.99|) in the pressure evolution equation 

OA, 



(3.94), we obtain 

3 dpr 1 d 



2 dt 



^m\H- Qt ' v ^ + ^ 2v0 dt 



Vdip 



(3.100) 



While there is no net turbulent heating, this does not rule out the possibi- 
lity of significant turbulent energy exchange between species (the turbulent heating 
species by species is the same order as the turbulent heat transport term, for in- 
stance). Few studies have been conducted investigating the effects of turbulent 
energy exchange [191 12Q] , and in these studies, the turbulent heating was defined as 
((53 ■ <5E)). Examining Eqs. (3.80)-(3.89), we see that this is not quite equivalent 



to the turbulent heating defined here (the difference between them is the additional 



term appearing in Eq. (3.85)). Consequently, the impact of the turbulent heating 



term in Eq. (3.93) on the evolution of equilibrium pressure profiles deserves further 
study. 



As an aside, we note that the net turbulent heating defined in Refs. [19] and [20] 
is also zero. From Poynting's Theorem (using low-frequency Maxwell's equations 
for the fluctuating fields), we have 



0_ 
dt 



5B< 



+ 



57T 



Air 



jdS ■ (5E x 5B) 



- / d 3 r 53 ■ 5E. 



(3.101) 



Using statistical periodicity of the fluctuations to eliminate the surface integral and 
applying the intermediate time average, we find 



cf r 53 ■ 5E = 0. 



(3.102) 
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Therefore, net turbulent heating of the equilibrium defined in this alternate way (as 
j d 3 r S3 ■ <5E) is also zero. 

In numerical simulations, it is often the case that the turbulent heating term 



as written in Eq. (3.93) has large amplitude oscillations in time, making it difficult 



to quickly obtain a steady-state time average. Before we finish our derivation, we 
would like to rewrite the turbulent heating term in a more convenient form for 
simulation. Making use of the fact that h t is independent of gyroangle at fixed R, 
we can change variables from r to R in our phase space integration to get 



d3r f J3 u d * 

~V qht m 



d 3 R 



d 3 v qh, 



dt ' 



(3.103) 



Multiplying the gyrokinetic equation (3.46) by h t T /F , averaging over phase space 
(with R as our spatial variable), and averaging over the intermediate time scale, we 
find that most terms do not contribute at this order. We are left with 



d v g-h t 



h t T 
F 



v x • VF - (C[h t ))^ 



(3.104) 



The first term on the righthand side in Eq. 3.104 is the energy exchange between 
the equilibrium and the turbulence, which is generally cooling the equilibrium. The 
second term is the collisional heating, or entropy generation, which is a positive- 
definite quantity. This sign-definiteness facilitates quick calculation of converged 
steady-state values for the turbulent heating. 



Using these results in Eq. (3.93), we obtain the final form of our equation for 



the equilibrium pressure evolution for each species: 
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(3.105) 



3.3 Summary 

In this chapter, we began with the Fokker-Planck equation and introduced a set 
of ordering assumptions that allowed us to derive a closed set of equations describing 
the self-consistent evolution of turbulence and thermodynamics processes (transport 
and heating). The ordering assumptions are a variant of the standard 5f gyrokinetic 
ordering. In particular, we assumed that all quantities could be split into (well- 
separated) slowly and rapidly varying parts in space and time, choosing a definite 
space-time ordering in terms of the expansion parameter p/L. The amplitude of 
of the rapidly varying parts were assumed much smaller than that of the slowly 
varying ones. Furthermore, we chose a definite ordering for the scale of velocity 



42 



space structures (5v ~ v t h) and for the collision frequency (v ~ uj). Our ordering 
and procedure for deriving the hierarchical equations follows closely the treatment 
of Refs. (13] and [IB]. A similar hierarchical set of equations, with the inclusion of 
low Mach-number, long- wavelength flows, is derived in Ref. [48J. 

For convenience in later chapters, we collect the key results from our calcula- 
tions here. The single-particle distribution function is written in the form 



f = F + h + f 2 , 



(3.106) 



where h is the non-Boltzmann part of the lowest-order perturbed distribution func- 
tion, /2 represents higher-order corrections, and 



F = F M (R) exp 



(3.107) 



with Fm a Maxwellian in velocity space. The slowly-varying, equilibrium part of the 
distribution function, F M , is independent of gyroangle and does not vary spatially 
within a magnetic flux surface: 



OF, 



b • VF M 



0. 



0. 



(3.108) 
(3.109) 



The non-Boltzmann part of the lowest-order perturbed distribution function (rep- 
resenting rapid fluctuations in space and time), h, is also independent of gyroangle 
at fixed guiding center position, R: 



dh\ 



- 



o. 



(3.110) 



The evolution of the part of h associated with turbulent dynamics is given by the 
well-known gyrokinetic equation: 

^ +vfi ■ Vh + (v x ) R • V(F + h) + v B • Vh = (C[h t ]) R + ^ d -^- (3.111) 

To close the system, we need information about the electromagnetic fields and the 
evolution of equilibrium (thermodynamic) quantities. In the gyrokinetic ordering 
used here, the low-frequency Maxwell's equations become: 
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(3.114) 
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The first expression above is Poisson's equation, with the assumption of quasineu- 
trality, and the remaining expessions are the parallel and perpendicular components 
of Ampere's Law. Finally, the equations describing the evolution of the equilibrium 
thermodynamic quantities (assuming no equilibrium flows) are 



dn 

dt ~ V dip \ 



1 } //V I rf 3 v,R 2 V0- 



—^F + h t V X C[h, 

ot q 



(3.115) 
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Chapter 4 

Resolving velocity space dynamics in gyrokinetics 



4.1 Introduction 

Many plasmas of interest to the astrophysical and fusion communities are 
weakly collisional. For such plasmas, velocity space dynamics are often impor- 
tant, and a kinetic description is necessary. Since the kinetic description requires a 
six-dimensional phase space, simulating weakly collisional plasma processes can be 
computationally challenging. Employing the gyrokinetic ordering [HI El EE]] reduces 
the dimensionality by eliminating gyrophase dependence, but we are still left with 
a relatively high-dimensional system. Consequently, one would like to know how 
many grid points are necessary along each dimension, particularly in velocity space, 
in order to resolve a given simulation. 

In the absence of collisions or some other form of dissipation, the distribution of 
particles in velocity space can develop arbitrarily small-scale structrues. Clearly, this 
presents a problem for gyrokinetic simulations, as an arbitrarily large number of grid 
points would be necessary to resolve such a system. Of course, all physical systems 
possess a finite collisionality, which sets a lower bound on the size of velocity space 
structures and, therefore, an upper bound on the number of grid points required 
for resolution. We would like to know how sensitive the plasma dynamics are to 
the magnitude and form of the velocity space dissipation. In particular, we would 
like answers to the following set of questions: Given a fixed number of grid points, 
how much dissipation is necessary to ensure a resolved simulation? Alternatively, 
given a fixed amount of dissipation, how many grid points are necessary to ensure a 
resolved simulation? Futhermore, what measurable effect, if any, does the addition 
of dissipation have on collisionless plasma dynamics? 

These questions have been addressed for very few plasma processes |l9j [50] , 
in large part due to the computational expense involved with such a study. In this 
chapter, we propose computationally efficient diagnostics for monitoring velocity 
space resolution and apply these diagnostics to a range of weakly-collisional plasma 
processes using the continuum gyrokinetic code GS2 [51]. With the aid of these 
diagnostics, we have implemented an adaptive collision frequency that allows us to 
resolve velocity space dynamics with the approximate minimal necessary physical 
dissipation [52J. We find that the velocity space dynamics for growing modes are well 
resolved with few velocity space grid points, even in the collisionless limit. Including 
a small amount of collisions [y <C uj) is necessary and often sufficient to adequately 
resolve nonlinear dynamics and the long-time behavior of linearly damped modes. 



The chapter is organized as follows. In Sec. 4.2 we discuss velocity space 



dynamics in gyrokinetics and provide examples illustrating the development of small- 
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scale structure in collisionless plasmas. Sec. 4J3 contains a brief overview of the 
Trinity velocity space grid and its dissipation mechanisms. We describe diagnostics 
for monitoring velocity space resolution in Sec. |4.4 and apply them to a number of 
plasma processes. In Sec. 4.5, we introduce an adaptive collision frequency and 



present numerical results. We discuss our findings in Sec. 4.6 



4.2 Gyrokinetic velocity space dynamics 

Trinity solves the coupled system consisting of the low-frequency Maxwell's 
equations and the nonlinear, electromagnetic gyrokinetic equation with a model 
Fokker-Planck collision operator: 



Oh 
dt 



K 



+ hib + v x + v B j • vh = (c[h]) n + T 
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where 



h 



is the non-Boltzmann part of the perturbed distribution function, 



n 



X 



b- V b + 



2 B 



is the sum of the curvature and VB drift velocities, 

c 



x 



b x Vx 



o 



is the generalized E x B velocity, 



X 



A 



(4.1) 



(4.2) 



(4.3) 



(4.4) 



(4.5) 



is the generalized electromagnetic potential, (-) R denotes a gyro-average at fixed 
guiding center position R, and 



M 



1 - 



(4.6) 



is the lowest order expansion of a Maxwell-Boltzmann distribution. The exact form 



of the collision operator, C[h], used in Trinity is discussed briefly in Sec. 4.3 and 
described in detail in Chapters 5 and 6. 



We can group the various terms in the gyrokinetic equation (4.1) into three 



distinct categories: source terms, labeled by S, which typically drive large-scale 
structures in velocity space; convection terms, labeled by K,, which lead to phase- 
mixing and the development of small-scale structures in velocity space; and dis- 
sipation, given by the collision operator, which smooths the distribution function 
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towards a shifted Maxwellian velocity distribution. In general, the structure that 
develops from the balancing of these terms can be quite complicated. However, 
we can gain insight into how small-scale velocity structures develop by considering 
simplified collisionless systems. 

In the absence of collisions, arbitrarily small scales can develop in velocity 
space. This is a result of phase-mixing, arising due to convection in real space [53j 
[15] . As a simple example of this phenomenon, we include in Appendix [B] a calcu- 
lation of the perturbed distribution function for the collisionless ion acoustic wave 
in a slab. The result, quoted here, illustrates the tendency of collisionless plasma 
processes to drive small-scale velocity space structures: 

fx (z, v }l , t) = ) + H(z, „„ , t) , (4.7) 

where the overbar on f\ indicates an average over perpendicular velocities. The 
quantities G and H are explicitly derived in Appendix [B] Here, it is sufficient to 
note that both G and H are smooth functions of the parallel velocity. The presence 
of the oscillatory factor e~ lk ^ v ^ 1 in the first term (often called the ballistic term) leads 
to the development of a characteristic wavelength in velocity space that decreases 
inversely with time. The amplitude of this ballistic term remains comparable to the 



second term in Eq. (4.7) for all time, leading to the development of large amplitude 
oscillations of the distribution function at arbitrarily small-scale s in velocity space. 
A snapshot of this behavior at t = 10 (k\\v t! i) is shown in Fig. 



4.1 



The same calculation carried out for the collisionless ITG mode in a slab yields 
a distribution function with a similar ballistic term component. However, since this 
mode is linearly unstable, there is also a term describing large-scale structure in 
velocity space whose amplitude grows in time to dominate the distribution function. 
As a result, no significant small-scale structure develops. This is a typical feature 
of linearly growing modes in the collisionless limit. 

Of course, all physical systems have a finite collisionality. The dissipation 
arising from this collisionality is critically important: It is a necessary requirement 
for the existence of a statistically steady state [5U [53], and it sets a lower bound 
on the scale-size of structures in velocity space [15]. A simple estimate for the 
scale-size of velocity space structures can be obtained by assuming a steady state 
and balancing the collisional term with the other terms in the gyrokinetic equation. 
Noting that C ~ vv\ h d 2 v (see e.g. Ref [21] or Ref [25]), we find 

Sv 117 , 

J- 4.8 

Vth V ^ 

where v is the collision frequency, u is the dynamic frequency of interest, Vth = 
^2T/m is the thermal velocity, and 5v is the scale-size of fluctuations in velocity 
space. This estimate predicts that velocity space structures much smaller than the 
thermal velocity develop in the weakly collisional limit, v <C u, as we would expect 
from our consideration of simplified collisionless systems. 
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Ion acoustic wave distribution function 




v 



Figure 4.1: Plot of f(v») (normalized by Fo) at t — 10 (k\\Vt,i) ■ The parallel veloc- 
ity on the horizontal axis is normalized by Vth and f(v\\) was initially a Maxwellian. 
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4.3 Trinity velocity space 



In order to understand the velocity space resolution diagnostics described in 
later sections, it is necessary for the reader to have a basic knowledge of the way 
in which velocity space dynamics are treated in Trinity. To that purpose, we now 
give a brief explanation of the velocity space coordinates and dissipation mechanisms 
employed in Trinity. 

4.3.1 Velocity space coordinates 

Only two velocity space coordinates are necessary in gyrokinetics because gy- 
roaveraging has eliminated any gyrophase dependence. Fundamentally, Trinity 
uses kinetic energy, e, and a quantity related to magnetic moment, A = fi/e, as its 
velocity space coordinates. This choice eliminates all velocity space derivatives from 
the collisionless gyrokinetic equation and simplifies the discretization of derivatives 
in the model collision operator. Consequently, the spacing of the velocity space 
grid points is chosen to provide accurate velocity space integrals while satisfying the 
necessary boundary condition at particle bounce points. 

4.3.1.1 Energy grid 

The volume element in velocity space can be written 

where •& is the gyroangle and a denotes the sign of v». Until recently, the energy 
grid in GS2 followed the treatment of Ref. [55], which places energy integrals in a 
convenient form by a change of variables to 

X(x) = — %xe~ x2 + Erf [x] , (4.10) 
v 71 

where x = v/vth- This transforms the range of integration from x G [0, oo) to 
■V ... 0.1): 

/r~ f2ir pi/Bo j\ pi 

^=f*4E/ v^ml dXe "- (4U) 

The integration domain is split into the subintervals [0, Xq) and [X ,l), with the 
perturbed distribution function assumed to be approximately Maxwellian on [Xq, 1). 
Gauss-Legendre quadrature rules [56] are then used to determine the location of the 
grid points in the interval [0,Xq). 

This energy grid provides spectrally accurate energy integrals (i.e. error ~ 
(1/N) N , where iV is the number of energy grid points), provided the integrand 
is analytic over the integration domain (see e.g. Ref. |57J). Unfortunately, this 
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Figure 4.2: Plot of normalized velocity x over the entire X domain. The func- 
tion x(X) has singularities at the boundaries of the domain due to a branch cut 
originating at X = and to x going to oo at X = 1. 



is seldom the case. To understand why, we consider the functional form of x(X). 
Taylor expanding X about x = 0, we find X ~ x 3 , or equivalently, x ~ X 1 / 3 . This 
indicates a branch cut in x originating from X = 0, so that most functions of x are 
non- analytic at X = 0. In Fig 4.2 we examine x(X). We see that not only is x 
non-analytic at X = 0, but also at X = 1, where x — > oo. The fact that x possesses 
singularities at the endpoints of the domain in X means that the integration scheme 
is not spectrally accurate for most integrands of interest (especially since the Bessel 
functions Jo(k±v±/Q) and Ji(k±v±/Q), which are non-analytic at X = and X = 1, 
appear in all integrals of the distribution function at fixed particle position r). This 



is demonstrated in Fig. 4.3 where we examine the accuracy of the numerical integral 
of h(R) = Fm (at fixed r) as we vary the number of velocity space grid points. 

In order to achieve spectral accuracy, we have implemented a new energy grid. 
We begin by splitting the velocity integration into two separate integrals: 



dx x 2 G(x) 



dx x 2 G(x) + dx x 2 G(x 



(4.12) 



.I'd 



where x$ is a free parameter and G(x) is the function we wish to integrate. On the 
first interval, (0,xo), we use Gauss-Legendre quadrature rules in x to obtain grid 
locations. Note that use of x as our integration variable ensures that the integrand 
x 2 G(x) will be analytic as long as G is analytic in x over the interval. 
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Figure 4.3: Plot showing absolute error in numerical integral of Jq{ q" x )- The inte- 
gration scheme of Ref. [55] has error proportional to negrid~ 3 ' 25 , while our scheme 
has error proportional to 0.6 * negrid~ os * ne9nd . Note that the minimum error in the 
non-spectral scheme is on the order of 10~ 6 , while in our scheme it is on the order 
of 10~ 16 , which is a limitation imposed by double precision evaluation of the Bessel 
function. 
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For the interval (xq, oo) we make the change of variable y 
the integral to 



X Xt 



dx x 2 G(x) 



dy e 



! \Jy + xlG(x) 



to transform 



(4.13) 



We then use Gauss-Laguerre quadrature rules in y to obtain grid locations. Note 
that the volume element is analytic within the domain of integration, as is x(y) = 
a/xq + y, so that the integrand will be analytic as long as G is an analytic function 
of x. 

Our use of spectral integration techniques (i.e. Gaussian quadrature), coupled 
with the analyticity of our integrand for well-behaved functions G(x), ensures the 
spectral accuracy of our integration scheme. While an exponential order of con- 
vergence is assured, the rate of convergence depends on the exact nature of the 
integrand and our choice of the parameter xq. In general we choose xq > 2.5 so that 



the branch cut at y 



.X c 



is sufficiently far from the domain of integration in y to 



minimally impact the rate of convergence. We demonstrate the spectral accuracy of 
the scheme and determine the rate of convergence for a number of test functions in 



Figs. 4.3 and 4.4 It is worthwhile to note that for few grid points (< 8 in Fig. 4.3) 
the grid given in Ref. [55] may be more accurate. This is because the energy variable 
X eliminates velocity-dependence of the volume element (when solving for the nor- 
malized distribution function h = h/Fo), while the new v-space integrals described 
here have the velocity-dependent volume element x 2 e~ x that must be integrated 
regardless of the form of h. 



4.3.1.2 Lambda grid 

For systems with curved magnetic field lines, special care is also required 
when dealing with A [51] . There are two reasons for this: the grid points provided 
by Gaussian quadrature rules are concentrated near the endpoints of the domain, 
whereas one would like them to be concentrated at the trapped-passing boundary; 
and one must ensure that the proper boundary condition (i.e. f(vu = + ) = f(vu = 
- )) is satisfied at each of the bounce points. Consequently, the A-grid is divided 
into two regions corresponding to trapped and untrapped particles, respectively. 

For values of A such that < A < 1/B max , the corresponding particles are 
untrapped by the magnetic potential well. In this region of velocity space, the 
integration variable £ = y/\ — \B max is chosen. It is similar to pitch-angle, but 
it has no spatial dependence. Similarly to the energy, Gauss-Legendre quadrature 
rules are used to obtain the location of grid points in £. This naturally provides a 
concentration of gridpoints near the trapped-passing boundary. 

For values of A such that 1/B max < A < 1/B min , the corresponding particles 
are trapped by the magnetic potential well. In the trapped region, grid points are 
chosen to fall on bounce points in order to allow for the enforcement of boundary 
conditions. Mathematically, this means that for each value of 6, there must be a 
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Error in numerical integration of test functions g(v) 




-| e _-|4 I i i i i i i i i I 
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negrid 



Figure 4.4: Plot showing absolute error in numerical integral of a number of test 
functions, g(v) . Our integration scheme has a rate of convergence proportional to 
approximately O.Q*negrid~°' 6 * ne9nd . Note that the minimum error approaches 10~ 16 , 
which is a limitation imposed by double precision arithmetic. 
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GS2 velocity space grid 
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Figure 4.5: Typical velocity space grid used in Trinity. Grid points are concen- 
trated near the trapped-passing boundary (whose location varies with 9) and at 
lower energy values where the Maxwellian weighting dominates. 



corresponding A such that 



m(9) 



f (0) = J^-L = y/l - \B (6) = 0, (4.14) 

where 9 gives the position along the unperturbed magnetic field line and £ is the 
pitch-angle. This choice of A values also leads to a concentration of grid points near 
the trapped-passing boundary. A typical Trinity grid layout for a system with 



trapped particles is shown in Fig. |4.5| It should be noted that the A integrals, like 
the energy integrals, are spectrally accurate, provided the distribution function is 
analytic in £. 

4.3.2 Velocity space dissipation 

Some form of dissipation is often necessary to prevent the formation of arbi- 
trarily small-scale structures in velocity space. This can be achieved either through 
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artificial numerical dissipation or through implementation of a model collision op- 
erator. Both options are available in Trinity. 



4.3.2.1 Model collision operator 

Trinity uses a model Fokker-Planck collision operator that includes the ef- 
fects of pitch-angle scattering and energy diffusion while satsifying Boltzmann's 
H-Theorem and conserving particle number, momentum, and energy [2"T| |2"2"]: 

C[h]=£[h]+V[h]+M[h], (4.15) 

where 

„r, n vd ( d , ~x dh 1 d 2 h\ 

£ N = T SfC-^K + iiFw (416) 



is the Lorentz collision operator, 



p i ft i = hi {^ F 'mk) (417) 

is the energy diffusion operator, and M. [h] contains momentum- and energy-conserving 
corrections. We defer a detailed discussion of the collision operator to Chapter 5. 
Here we simply present the gyroaveraged collision operator in spectral form: 

~ Tea ~ Vd + € ) + p s{l-£ ) I h k + v E v J {a)F — — 

8fig Wfc / J d 3 v u E v 4 F 

n ( r i \ I d3v ^DV^hia)^ Jd 3 v v D v±Ji(a)h^ 

V J d 6 v v D v»Fo J d 3 v udv^F 

( f d 3 v Auv\\J (a)h k f d 3 v AuvxJ^ajh^ 

- AuF J (a )v\\ f " 2 y + J x {a)v x J \ J 

\ J d 6 v Auv^o J d 6 v Avv^Pq 

(4.18) 

where k is the wavenumber and a = kv±/Qo- The velocity-dependent collision 
frequencies u s , i>d, ve, and Av are given by 




(4.19) 

"d = ( ""^P ~ '-j ) (4.20) 



Av = v D - u 8 , (4.22) 

with v su the frequency of collisions of particles of species s with particles of species 
u. Details on numerical implementation of the collision operator (4.18) are given in 
Chapter 6. 
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4.3.2.2 Numerical dissipation 



Numerical dissipation enters in Trinity through two mechanisms. The first 
is the optional decentering of spatial and temporal finite differences, as described in 
Ref. [51]. The lowest order contribution to dissipation due to decentering in time 
and space is 



d 2 h 
dtd6 



AO 5 



h) 



i+1/2 



At [(3 



d 2 (X) n 
dtdO 



A9 5 



qFo 
T 



(4.23) 

where AO is the grid spacing along the field line, At is the time step size, and S and j3 
are parameters that allow for the variation of the spatial and temporal discretization 
schemes between fully explicit (5 or j3 = 0) and fully implicit (5 or j3 = 1)Q 

In order to see how this term leads to dissipation, we consider the simplified 
system governed by the equation 



dh dh 
di + V d6 



0. 



(4.24) 



Finite differencing this equation using the scheme given in Ref. [51] . we find that 
numerically we are solving the equation 



dh dh 
di +V d6 



d 2 h 
dtdO 



A9 5 



vAt p 



Assuming h = h(t)e l , we obtain the solution 



h(t) ~ exp 



kvt 



_i — k (AO (5 - 1/2) + vAt ((3 - l/2))_ 



(4.25) 



(4.26) 



which is damped unless (5 — 8 — 1/2, as show in Fig. 4.6[ While decentering of 
finite differences can sometimes improve numerical stability, care must be taken 
to ensure such artificial dissipation does not lead to unphysical behavior. This is 
typically done by monitoring the ratio of artificial to physical dissipation, which, 
ideally, should be small. 

The second source of numerical dissipation arises in systems with sheared 
magentic fields due to the necessity of a 'twist-and-shift' parallel boundary condi- 
tion [TT]. This non-periodic boundary condition couples modes at opposite ends of 
the simulation domain along the field line. Since only a finite number of modes can 
be kept in a simulation, some modes will eventually couple to modes that are not 
present, and this information is lost. The information that is lost is replaced by a 
smoothed distribution function, leading to a loss of entropy. If a sufficiently large 
number of modes are kept in the simulation, the energy contained in the highest 
modes, and therefore the lost entropy, should be negligible. Currently, this entropy 
loss is not regularly diagnosed in Trinity. In principle, it could (and should) be 
diagnosed in order to verify that the entropy lost is small compared to the entropy 
generated by collisions. 



1 Trinity actually uses (3 = (3 — 1/2, but we choose to use (3 here for simplicity. 
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Figure 4.6: Damping of the real part of the distribution function h [Eq. (4.26)] as 
a result of decentered finite differences in space and time. Here, we are considering 
v = 1, k = 2, Ax = At = 0.1, and (5 = 5 = 1.0 (fully implicit). 
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4.4 Velocity space resolution diagnostics 



There are numerous ways in which one could try to determine whether or 
not a particular simulation is well-resolved in velocity space. Ideally, one would 
perform a grid convergence study for each simulation; if quantities of interest are 
unchanged by doubling the number of grid points, one can feel relatively confident 
in the simulation results. However, this process is computationally expensive, as 
it involves running a simulation multiple times with an excessive number of grid 
points. Consequently, it is not desirable to perform a grid convergence study for 
every simulation. In practice, one tests convergence for a problem thought to be 
resolution intensive and posits that other simulations, which likely require fewer grid 
points, are therefore resolved. Unfortunately, one seldom knows in advance how fine 
the structure in velocity space will become, so one can't be fully confident that every 
simulation is resolved. 

An alternative approach that has recently gained popularity in the computa- 
tional plasma physics community involves monitoring entropy balance in the sys- 



tem [49, 50J. Multiplying the gyrokinetic equation (4.1) by hT /F and integrating 



over all phase space gives the desired relation for entropy balance: 

1 (99 

^ = * + r + ff, (4.27) 

where 

S = j d 3 r J d 3 v ^h 2 (4.28) 
is a lowest order expression for the perturbed entropy, 

X = J d 3 r J d\ q^^h (4.29) 

describes turbulent heating, 

r = - J d 3 r J d\ /iT v x • V In F (4.30) 
is the entropy flux due to background inhomogeneity, and 

J3„ f j3_ hT o 



H = J d 3 r J d 6 v (C[h]) (4.31) 

describes entropy change due to collisional heating. 

Since the gyrokinetic equation itself is automatically satisfied by a gyrokinetic 



solver, the only possible sources of inbalance in Eq. (4.27) come from numerical 
dissipation and errors in the numerical approximations to phase space integrals. If 
the change in entropy due to numerical dissipation is also diagnosed and included 
in the entropy balance, as is often the case, then we are left with errors due only 
to phase space integration. Since the errors in these particular integrals are not di- 
rectly related to errors in the calculation of the distribution function at the newest 
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timestep, they do not necessarily correlate with the simulation resolution. In par- 
ticular, one could easily define a poorly-resolved system for which this diagnostic 
predicts perfect entropy balance. One such example is the linear, collisionless ion 
acoustic wave in a slab (treated in detail in Appendix [Bj. For this case, we numer- 
ically find entropy balance despite the fact that the numerical damping rate goes 
bad due to poor resolution in velocity space. 

Of course, one could simply produce plots or movies of the distribution function 
in velocity space over the course of the simulation to see if structure develops at the 
gridscale. This is undoubtedly useful and possibly sufficient in some cases. However, 
what exactly one sees depends on how the data is visualized; for data on irregularly 
spaced grids, the interpolation scheme used to generate the images often introduces 
erroneous or misleading structure. Furthermore, for simulations involving non-trivial 
spatial structure, one would have to examine movies of the distribution function at 
each point in physical space. This is a memory- and time-intensive approach that 
is rarely feasible. 

We would like to have computationally cheap diagnostics that provide real- 
time information on velocity space resolution that is easy to analyze and interpret. 
In the following subsections, we present two such diagnostics developed for imple- 
mentation in Trinity that could easily be adapted for use in other continuum kinetic 
simulations. 



4.4.1 Integral error estimates 

Upon consideration of the collisionless gyrokinetic-Maxwell's system of equa- 
tions, one finds that the only nontrivial operation in velocity space is integration, 
which enters in the calculation of the electromagnetic fields. Consequently, resolu- 
tion in velocity space is limited only by the accuracy with which the velocity space 
integrals are calculated. By calculating the error in our numerical integration, we 
are thus able to monitor velocity space resolution. 

In particular, when we discretize the gyrokinetic equation, we obtain an equa- 
tion of the form 

9j+i = G \9j, $ i> Xj, Xj+i] , (4-32) 

where g = (fi) is the perturbed, guiding center distribution function, $ is the 
electrostatic potential, \ is the generalized electromagnetic potential defined in 



Eq. (4.5), G is a function that depends on the details of the numerical scheme, 
and the subscript denotes the timestep. We assume that the time-converged solu- 
tion for g is independent of the initial condition. Since using the calculated gj and $j 
is equivalent to specifying a new initial condition, we find that the time-converged 
solution is independent of errors in g and $ at earlier timesteps. This is convenient 
because it means we can monitor resolution merely by calculating the error made 
in the latest timesteps of a time-converged simulation. 

Ideally, we would accomplish this by calculating estimates for the error in 



and Xj+i an d plugging these into Eq. (4.32) to obtain an error estimate for 



gj+i- This might be feasible for linear systems, but the presence of nonlinear terms 
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makes this approach computationally prohibitive. Consequently, we must define an 
alternative quantity whose error estimate is cheaper to compute, but that can still 
be used as a means of monitoring velocity space resolution. There are numerous 
possible candidates; we choose to compute two quantities, t>$ and va, related to 
V±$ and V±Af 

where k x and k y are the wavenumbers corresponding to the x — y coordinates x = 
(ip — ipo) qo/B r and y = — (a — oq) r /q [11]. Here, ip is the poloidal flux, a is 
the field line label, B is the background magnetic field at the magnetic axis, r is 
the distance from the magnetic axis to the center of the simulation domain, and go 
is the safety factor on the field line of interest, labeled by (^o,«o)- The quantities 



in Eq. (4.33) were chosen because, with the exceptions of the parallel convection 
term and one source term, $ and An always enter the gyrokinetic equation for 
g multiplied by either k x or k y . Therefore, it is reasonable that this k- weighted 
quantity is most likely to be responsible for errors in gj+\. Although not considered 



here, the expression (4.33) could potentially be improved by including k\\ in the 
max operator. This would take into account the effect of the parallel convection 
term. However, there is recent theoretical [15] and numerical [HI [T7j evidence 
which suggests that velocity space structure is primarily generated by nonlinear 
perpendicular phase mixing (instead of linear, parallel phase mixing). 

Having chosen appropriate indicators of velocity space resolution, we must 
devise a method for estimating the error in these quantities. This error depends on 
the particular numerical integration scheme used. For the energy and untrapped A 
integrals, which use Gaussian quadrature, the error, e<3, is given by 

e G = 7«J (2m) (C), (4-34) 

where / is the integrand, m is the number of grid points, and ( is some unkown 
point in the interval of integration. The quantity 7 m is 

2 2m+1 (ml) 4 , . 

7m = — - — * 4.35 

(2m + l)[(2m)!] 3 V ' 

for the untrapped A and finite domain energy integrals that use Gauss-Legendre 
quadrature and 

(ml) 2 

for the semi-infinite domain energy integral that uses Gauss-Laguerre quadrature. 
The error, e^, for the trapped A integrals, which use a newly upgraded integration 
scheme based on Lagrange interpolating polynomials (see, e.g. Ref. [SS]), is given 
by 

£l = — } [ f (m \(Mx)dx, (4.37) 
ml / 
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where 



nix) = n 



(4.38) 



8=1 



with Xi the i th grid point. It should be noted that ( in Eqn (4.37) is an unknown 



function of x whose domain is some subset of the interval of integration. 

From Eqns (4.34) and (4.37), we see that Gaussian quadrature gives exact 



results for polynomials of degree less than 2m, while the Lagrangian method gives 
exact results only for polynomials of degree less than m. We say that the two 
schemes have degrees of precision 2m — 1 and m — 1, respectively. This difference 
arises because the grid points in the Lagrangian method are fixed by boundary 
conditions, whereas the grid points in Gaussian quadrature are free parameters 
optimally chosen to improve the scheme's degree of precision. 



Unfortunately, the formal error expressions (4.34) and (4.37) are not very 



useful in practice: they require information about high-order derivatives of the dis- 
tribution function, which is unavailable. As an alternative estimate for the error, we 
choose to compare multiple integral approximations computed with different degrees 
of precision, a common technique in numerical analysis. 



4.4.1.1 General description of the scheme 

Given the value of a function f(x) at N fixed points on the interval [a, b], 
we would like to find two different approximations to the integral J f(x)dx. In 
our earlier discussion, we stated that an approximation with degree of precision 
iV — 1 can be found using a technique based on Lagrange interpolation; we call this 
approximation Ah- If we instead choose to use only M of the given functional values 
(M < N), we can use the same technique to find another integral approximation, 
Ai, with degree of precision M — 1. An estimate for the absolute error e a in the less 
accurate of these two approximations is obtained by taking the difference between 
the two: 

e a =\A h -Ai\. (4.39) 

Making the reasonable assumption that the approximation with higher degree of 
precision is more accurate, e a represents the error in A\. However, it can also be 
used as a more conservative error estimate for Ah- 

If the N points are chosen according to Gaussian quadrature rules, then one 
can find an integral approximation with degree of precision 2N — 1. As before, 
a second approximation can be obtained by using only M of the N grid points. 
However, due to the uniqueness of the grid points used for Gaussian quadrature, 
the M-point grid no longer satisfies Gaussian quadrature rules. As a result, this 
second approximation once again has degree of precision M — 1. Since the degrees 
of precision of the two approximations differ by greater than a factor of two, the 
resulting error estimate is likely to be very conservative when applied to Ah- The 
factor of approximately two difference in degree of precision makes this error estimate 
similar to that obtained by comparing results from runs with N and N/2 grid points, 
respectively (for which the degrees of precision would be 2N — 1 and N — 1). 



61 



The conservative nature of the error estimate for Ah depends upon our assump- 
tion that a higher degree of precision results in a more accurate integral approxima- 
tion. For Gaussian quadrature, it can be shown that the error in the integral ap- 
proximation can be made arbitrarily small by choosing the degree of precision large 
enough [56J. The same result does not necessarily hold for the Lagrangian method 
with arbitrary grid spacing because the weights in this case are not all guaranteed 
to be positive. However, the error €m in an M-point integral approximation satisfies 

e M <2 
< 



where e can be chosen arbitrarily small for large enough M, and w\ m> is the weight 
corresponding to the i th grid point out of M. From this result, we see that as long 
as k is bounded when M — * oo, then €m — ► as M —>■ oo. This cannot be verified 
in advance, but one can gain confidence by checking a posteriori. In practice, we 
calculate k for the chosen M and subdivide the integration domain into subintervals 
with fewer points if k is larger than some reasonable value. 



M 



(M) 



1=1 

2eM max \w; 

i=l,M 

2eAf«(M), 



(4.40) 

(4.41) 
(4.42) 



4.4.1.2 Implementation in Trinity 

In Trinity, we must compute two-dimensional integrals over energy and A. 



As stated in Sec. |4.3[ each of these integrals is effectively separated into two by 
splitting the A integration into trapped and untrapped regions. Since the number of 
grid points in energy and both A regions can be varied independently of each other, 
we wish to monitor resolution in each of these three variables individually. This 
entails computing three separate integral error estimates: one for energy integrals, 
one for untrapped A integrals, and one for trapped A integrals. 

These integral error estimates are calculated using the technique described in 
the previous subsection. For energy and untrapped A integrals, Gaussian quadrature 
is used to obtain the two-dimensional integral approximation Ah- This approxima- 
tion has degree of precision 2N e — 1 for the energy integration and 2N U — 1 for the 
untrapped A integration, where N e and N u are the number of energy and untrapped 
A grid points, respectively. To obtain the second approximation, Ai, we fix the grid 
and weights for one variable and drop one grid point for the other variable, recom- 
puting the weights. As an example, we choose to drop an untrapped A grid point. 
The degree of precision for A\ is then 2N e — 1 for the energy integration and N u — 2 
for the untrapped A integration. Since there is nothing special about the particular 
grid point we drop, we repeat the process a total of N u times, each time dropping a 
different point and computing a different set of weights. The final error estimate is 
an average of these error estimates. 

For the trapped A integrals, Lagrangian quadrature is used to obtain Ah, which 
has degree of precision N t — 1. We obtain the approximation Ai by dropping two 
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Figure 4.7: Red grid points are sample trapped A grid points that are dropped when 
calculating integral approximation with lower degree of precision. 



points symmetrically about v\\ = 0, as shown in Fig. 4.7 We drop an additional point 
here because it provides a slightly more conservative error estimate and because 
maintaining the symmetry of the grid points provides better stability for the weights 
associated with the Lagrange interpolation scheme. As before, we repeat this process 
for each possible grid point pair and take the average of the individual error estimates 
to get the final error estimate. 

All modified grids and weights necessary for the integral error estimates are 
computed once at initialization and need not be computed again. The additional 
integrations necessary to obtain our error estimates are computationally cheap when 
compared to the expense of solving for the distribution function and fields at each 
time step. Furthermore, we do not need an error estimate at each time step, so 
the diagnostic can be used sparingly. Consequently, our error estimate comes at 
essentially no extra cost. 
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4.4.2 Spectral method 



An alternative method for testing v-space resolution is to expand the velocity 
space distribution function in an appropriate basis set and monitor the amplitude 
of the basis function coefficients. Whenever the highest mode number coefficients 
that can be accurately calculated in the simulation acquire appreciable amplitudes, 
we can no longer feel confident that the simulation is resolved. Since we choose 
our grid points according to Gauss-Legendre quadrature, it is convenient (and most 
accurate) to choose the Legendre polynomials as our basis functions. The coefficient 
of the m th Legendre polynomial in the expansion of h is given by 



2 

2m + 1 



2m + 1 f 1 , . , . . , 

J h(s)P m (s)ds (4.43) 

n e — 1 

£ wM Si )P m (si), (4.44) 
i=i 



where P m is the m Legendre polynomial, and {w;} are the weights associated with 



Gauss-Legendre quadrature. The integral approximation in Eq. (4.44) has degree 
of precision 2N — 1. Assuming h has a degree of at least m (otherwise c m — 0), our 
approximation for c m is only exact for m < N. 

There are various ways in which one could use these {c m } to estimate the error 
in velocity space resolution. We assume locality of interaction between the various 
modes so that we only have to monitor the amplitudes of the few highest modes. At 
each (6, k x , fc y )-point, we find the maximum amplitude of the three highest mode 
number spectral coefficients, Ch, ma x, and the maximum amplitude of all the spectral 
coefficients, c max . We then use the following normalized sum as a relative estimate 
for the error: 

= ^ ^ Ch,max{@ )k x ,ky) j ^ ^ c max {6 , k Xl ky) . (4.45) 

When the normalized amplitude e c grows too large, we can no longer be confident 
that the simulation is resolved. Of course, how large e c can get before resolution 
suffers varies from problem to problem. As before with the integral method, we 
determine a scaled estimate of the error based on empirical evidence from a wide 
range of simulation data. 



4.4.3 Application of error diagnostics 

We have applied both the integral and spectral error diagnostics to a diverse 
set of simulations, including: linearly growing modes such as the electron drift 
wave and the ITG mode; linearly damped modes such as the ion acoustic wave and 
kinetic Alfven wave; neoclassical transport; and nonlinear dynamics of slab ETG 
and toroidal ITG modes. From these simulations, we have determined empirical 
scaling factors for our conservative error estimates. Here, we present typical results 
from a cross-section of the above simulations. 
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Figure 4.8: Comparison of actual and (unsealed) estimated error in wave frequency 
due to insufficient resolution in energy (top left), untrapped A (top right), and 
trapped A (bottom). The actual wave frequency, u, is determined from a higher 
resolution run with 64 grid points in energy and both trapped and untrapped A. 

The actual relative error, e, is then defined to be e = \/ ^7^ , where u n is the 

V M 

approximation to to obtained from a run with n grid points. 
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Figure 4.9: Non-Boltzmann part of the perturbed distribution function, normalized 
by Fo(a/pi). The use of a polar grid in velocity space minimizes the number of grid 
points necessary for resolution. 



Fig. 4J3 compares the unsealed error estimates in energy and A with the actual 
errors in growth rate as we vary the number of grid points in a linear simulation of 
the collisionless toroidal ITG mode (using Cyclone base case parameters [30]). The 
simulation remains well-resolved down to very few grid points, and, qualitatively, 
the error estimates agree well with the actual error. The error due to resolution 
in untrapped A is still small for as little as four grid points due to our choice of 
velocity variables, as illustrated by the snapshot of the distribution function shown 



in Fig. 4.9 



Figs. 4.10 and 4.11 show the damping of A\\ and the corresponding scaled error 
estimates for the simulation of a collisionless kinetic Alfven wave. The collisionless 



damping rate in Fig. |4.10| agrees with theory until sub-gridscale structure devel- 
ops in velocity space, at which point damping ceases. The onset of sub-gridscale 
structure corresponds to the peak in scaled error in Fig. |4.11| The addition of a 



small collisionality prevents sub-gridscale structure, as shown in Fig. 4.10 , where the 
damping rate of A\\ agrees well with theory indefinitely. This is accurately predicted 
by the error estimates of Fig. 4.12[ which never reach appreciable magnitude. 
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Figure 4.10: Barnes damping of the kinetic Alfven wave. In the absence of collisions 
(left), sub-grid scale structures develop in velocity space, and the damping rate 
goes bad. A small collisionality [y <C 7) prevents the development of sub-grid 
scale structures in velocity space, and the damping rate remains correct indefinitely 
(right). 
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Figure 4.11: Integral (left) and spectral (right) error estimates for the collisionless 
kinetic Alfven wave. 
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Figure 4.12: Integral and spectral error estimates correctly indicate that the weakly 
collisional kinetic Alfven wave simulation is well resolved. 



4.5 Adaptive collision frequency 

As stated earlier, we would like to know what combination of dissipation and 
grid spacing is necessary for a resolved simulation. One way to approach this prob- 
lem is to fix the dissipation and vary the number of grid points to find how many 
are required to get an accurate result. This is the general idea behind the error 
estimation diagnostics described in the previous section. However, if we wanted to 
use this approach to ensure that the simulation remains resolved, we would have 
to implement an adaptive grid, which is difficult to do for massive, multi-processor 
calculations. 

Instead, we choose an alternative approach: we fix the number of grid points 
and vary the dissipation until we have a well- resolved result. In particular, we 
have implemented an adaptive collision frequency in Trinity that allows for the 
independent variation of the collisionality associated with pitch-angle scattering and 
energy diffusion. Given an acceptable error tolerance for velocity space calculations, 
a scaled version of the integral error estimate described in the previous section is 
used to determine whether or not the simulation is well-resolved. The collision 
frequency is then adjusted using a feedback process until the scaled estimate of the 
error converges to within some pre-specified window of the desired error tolerance. 
In this way, the approximate minimum possible dissipation is used to achieve an 
acceptable degree of resolution in velocity space. 

Of course, the amount of dissipation necessary to resolve a simulation at a fixed 
number of grid points may be quite large if a coarse grid is used. Consequently, the 
collisionless dynamics may be modified. As a result, it is necessary to compare the 
converged collision frequency with dynamic frequencies of interest in the problem. 
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Figure 4.13: (Left): Normalized electron heat flux vs. time for a nonlinear simulation 
of ETG turbulence. Scaled estimates of the error in energy and A resolution increase 
during nonlinear saturation, but are kept within the specified error tolerance of 
0.01 with the use of an adaptive collision frequency. (Right): Collision frequency 
(normalized by k\\Vth,e) vs. time. 



As an example we consider a nonlinear simulation of electron temperature 
gradient (ETG) turbulence in slab geometry (i.e. straight background magnetic 
field). In the nonlinear phase, small scales are expected to develop in velocity 
space, potentially challenging numerical resolution. In Fig. 4.13 , we see that this is 
indeed the case. Our velocity space resolution diagnostics indicate that the errors in 
velocity space begin to increase sharply during the transition from linear instability 
to turbulence. However, our use of an adaptive collision frequency prevents the 
estimated error from exceeding the user-defined relative error tolerance (in this case, 
0.01). We see that the error remains on the threshold of the error tolerance, while 
the collision frequency for energy diffusion increases to a steady-state value of v m 
0.27 k\\v t h,ei which is well below the dynamic frequency in the system. Consequently, 
the collisionless dynamics are unaltered. 



4.6 Summary 

In this chapter, we discussed the development of small-scale structure in ve- 
locity space, presented a set of velocity space resolution diagnostics for use in gy- 
rokinetic simulations, and introduced an adaptive collisionality that allows us to 
resolve simulations with an approximate minimal necessary dissipation for a fixed 



number of grid points in velocity space. In Sec. 4.2 we demonstrated the tendency of 
collisionless plasmas to develop increasingly fine scales in the distribution of particle 
velocities and discussed the phase mixing processes that lead to such behavior. 
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In Sec. 4J3 we described the treatment of velocity space in the gyrokinetic 
code Trinity. We gave details on the choice of velocity space variables (energy 
and pitch-angle) and discretization scheme, which is chosen to minimize the error of 
the numerical integrals necessary to obtain the electromagnetic fields. This included 
presentation of a newly implemented energy grid, which provides spectrally accurate 
integrals over particle energies. Additionally, we gave a brief discussion of both the 
physical and numerical dissipation mechanisms available for use in Trinity. 

We discussed common approaches to monitoring velocity space resolution in 
Sec. 4A and the difficulties associated with each. We then proposed two new mea- 
sures of velocity space resolution and detailed implementation in Trinity. One of 
the proposed resolution diagnostics involves obtaining estimates for the error in field 
integrals by comparing numerical integrals obtained using integration schemes with 
differing degrees of precision. The other resolution diagnostic involves decomposing 
the perturbed distribution function into spectral components in velocity space and 
monitoring the amplitude of the spectral coefficients. Both diagnostics should be 
quite conservative. 

We then applied our resolution diagnostics to a number of example problems, 
including Landau damping of the ion acoustic wave, Barnes damping of the ki- 
netic Alfven wave, and linear instability of the toroidal ITG mode. We found that 
both diagnostics do well in qualitatively estimating errors due to limited velocity 
space resolution. Due to their conservative nature, an empirical scaling factor was 
necessary to obtain correct quantitative predictions. 



In Sec. |4.5| we coupled the error estimates from our resolution diagnostics 
with a model physical collision operator to develop an adaptive collision frequency. 
This adaptive collision frequency allowed us to resolve velocity space while using 
an approximate minimal necessary amount of dissipation. When using the adaptive 
collision frequency, one must monitor the ratio of the collision frequency to the 
dynamic frequency to ensure that one is still within the weakly collisional regime. 

In conclusion, we found that dissipation was not necessary to resolve linear 
instabilities, but it was necessary to resolve nonlinear dynamics and linearly damped 
waves. For the nonlinear cases considered here (slab ETG and toroidal ITG), the 
required collisionality for resolution obtained with the adaptive collision frequency 
was found to be no larger than the physical collisionality used in modern fusion 
experiments. 
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Chapter 5 

Linearized model Fokker-Planck collision operator for 
gyrokinetics: theory 



5.1 Introduction 

It has long been known that in many turbulent systems the difference between 
vanishingly small dissipation and no dissipation is striking, and that this can be 
linked theoretically to the non-interchangeability of limits t — > oo and v — > 0, where 
v is e.g. viscosity, resistivity or collision frequency. Turbulence transfers energy 
from scales at which it is injected into the system to scales where it is dissipated, 
leading to heating. When the dissipation coefficients are small, the system has to 
generate very fine-scale fluctuations in order to transfer the energy to scales at which 
dissipation becomes efficient. However, with finite u, there will always exist a scale 
at which the injected energy is dissipated. 

In plasma turbulence, all dissipation (meaning any effect that leads to irre- 
versible heating) is ultimately collisional, so the transfer of energy generally occurs 
in phase space — i.e., both in the position and velocity space (see extended dis- 
cussion of energy cascade in plasma turbulence in Ref. [H] and references therein). 
There are a number of specific mechanisms, both linear and nonlinear, that give rise 
to phase-space mixing E31 O EHl EOl SSI US]- It is the resulting large gradients 
in the velocity space that eventually bring collisions into play however small the col- 
lision frequency (such small-scale velocity-space structure has, e.g., been found and 
explicitly measured in gyrokinetic simulations [HI EU, |22l [16] ) . Thus, in any plasma 
turbulence simulation, some effective collisionality should be present to smooth the 
small-scale structure in the velocity. 

While one may take the view that the numerical grid can play the role of effec- 
tive collisions [50], we consider it to be a safer course of action to model collisional 
physics in a controlled fashion. In order to explain why, we would like to emphasize 
that, besides velocity-space smoothing, there is another key reason why collisions 
must be included. Collisions, through the dissipation of small-scale fluctuations 
in phase space, provide the link between irreversible plasma heating (macroscopic 
transport) and turbulence, so they are necessary in order for the system to converge 
to a statistically steady state. We shall now explain this statement. 

Consider the 5f kinetics detailed in Chapter 2. This model assumes that 
it is physically reasonable to split the distribution function into a slowly (both 
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spatially and temporally) varying equilibrium part and a rapidly varying fluctuating 
part: / = F + 5f. We saw in Chapter 2 that F is a Maxwellian distribution, 
F = (n /vr 3 / 2 u^ h ) exp(— v 2 /v 2 h ), where n is density, v th = (2T /m) 1 / 2 is the thermal 
speed, Tq is temperature and m is particle mass. This will be the case if collisions 
are not extremely weak (for the weakly collisional formulation of 5f gyrokinetics, 
see Ref. [ID]). One can show that the fundamental energy balance governing the 
evolution of the turbulent fluctuations is [SH HH S21 US UHl [621 [63l El] 




T 0s Sf ( 



(5.1) 

C[Sf s ]dvdr, 

where s is the species index, SS — — ff drdv5f 2 /2F is the entropy of the fluctu- 
ations, U = J dr (E 2 + B 2 )/87r is the energy of the (fluctuating) electromagnetic 
field, P is the input power (energy source of the turbulence), and C[Sf) is the lin- 
earized collision operator. In many types of plasma turbulence studied in fusion 
contexts, the input power P is proportional to the heat flux and it is the parameter 
dependence of the mean value of this quantity in the statistically stationary state 
that is sought as the principal outcome of the simulations. We can see immediately 
from the above equation that collisions are required to achieve such a steady state 
(as has been shown in numerical simulations [SH [53j HHl |65]) and that in this steady 
state, P must be balanced on the average by the collisional dissipation term. 

The key property of the collision operator required for this transfer of energy 
from turbulence to the equilibrium distribution to work correctly and, therefore, for 
the heat fluxes to converge to correct steady-state values, is that the collision term 



in Eq. (5.1) must be negative-definite: 

JJ 6 -^C[5f]drdv<0. (5.2) 

This ensures that heating is irreversible and that collisions cannot decrease entropy, 
the latter being the statement of Boltzmann's if -theorem [66]. Any spurious sink of 
entropy will adversely affect the balance between turbulent fluxes and dissipation, 
so it is clear that any model for collisional dissipation must respect the if-theorem. 

In view of the above discussion, we can formulate a reasonably restrictive set of 
criteria for any model collision operator: providing dissipation at small scales; obey- 



ing the if -theorem [Eq. (5.2)]; locally conserving particle number, momentum, and 
energy; and vanishing on a (local, perturbed) Maxwellian distribution. While these 
properties are analytically convenient, for numerical simulations the operator should 
also be efficiently implementable and carry these properties (at least approximately) 
over to the numerical scheme. 

The effect of small angle Coulomb collisions on an arbirtrary distribution func- 
tion was originally calculated by Landau [23J. In the 5f kinetic model we naturally 
consider the linearized Landau operator [HZ]. However, it is sufficiently complex that 
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it would exceed the limits on numerical resources that can be realistically expended 
on modeling the collisional physics. Consequently, several simplified model collision 
operators have been developed, both for analytical and computational convenience, 
that try to capture the qualitative essence, if not the quantitative detail, of the 
physics involved j68j EH [25]. This course of action is, indeed, eminently sensible: 
from Eq. (5.1 ), it seems plausible that, at least as far as calculating integral charac- 
teristics such as the turbulent fluxes is concerned, neither the exact functional form 
of the collision operator (provided it satisfies the criteria discussed above) nor the 
exact value of the collision frequency (provided it is sufficiently small) should be 
important. All we need is a physically reasonable dissipation mechanism. 

For these purposes, it has often been deemed sufficient to use the pitch-angle- 
scattering (Lorentz) operator, sometimes adjusted for momentum conservation [681 
IB"?] . However, in kinetic turbulence, there is no reason that small-scale velocity- 
space structure should be restricted to pitch angles. In fact, standard phase-mixing 
mechanisms applied to gyrokinetics produce structure in v\\ j59j 09], and there is 
also nonlinear gyrokinetic phase mixing that gives rise to structure in v±, which 
may be an even faster and more efficient process [601 HSl E]. Thus, a priori one 
expects to see small scales both in pitch angle and in the energy variable (£ and v). 
Indeed, it has been confirmed in simulations [22J that with only Lorentz scattering, 
structure rapidly forms at the grid scale in energy. Thus, a numerically suitable 
model collision operator must include energy diffusion. 

In this chapterQwe propose such an operator (other operators including energy 
diffusion have been previously suggested [211 [25] ; we include a detailed comparison 
of our operator with these in Appendix [E]). Our model operator for like-particle col- 
lisions, including both pitch-angle scattering and energy diffusion and satisfying all 
of the physical constraints discussed above, is given in Sec. 5.2| (the proof of the H- 
theorem for it is presented in AppendixjC]). In Sec. 5.3, it is converted (gyroaveraged) 
into the form suitable for use in gyrokinetic simulations — a procedure that pro- 
duces some nontrivial modifications. In Sec. 5.4, we explain how interspecies (and, 
in particular, electron-ion) collisions can be modeled in gyrokinetic simulations to 
ensure that such effects as resistivity are correctly captured. Section [5]5 contains a 
short summary and a discussion of the consequences of the work presented here. 

The anlytical developments presented in this chapter form the basis for the 
numerical implementation of collisions in the publicly available gyrokinetic code 
Trinity. This numerical implementation, as well as a suite of numerical tests are 
presented in Chapter 6. 



5.2 A New Model Collision Operator 

In this section, we present a new model collision operator for like-particle 
collisions that satisfies the criteria stated above. The interpecies collisions will be 
considered in Sec. 15. 41 

Let us start by introducing some standard notation. In discussing collision 

lr This chapter is taken from a co-authored paper currently in press |21j 
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operators on phase space, we shall denote r the position variable in the physical 
space and use the (v,£, $) coordinates in velocity space, where v = \v\ is the en- 
ergy variable, £ = v\\/v is the pitch-angle variable, and i? the gyroangle about the 
equilibrium magnetic field. One can easily adapt the operators presented here to 
unmagnetized plasmas, but as we are interested in gyrokinetic plasmas, we shall 
concentrate on the strongly magnetized case. Taking the notation of Ref. (67] as 
the standard, we introduce the normalized velocity variable x = f/i>th and a set of 
velocity-dependent collision frequencies for like-particle collisions: 

Erf(x) - G(x) 

v D {v) = v , 5.3 

x 6 

v a {v) = v — — , 5.4 

X 

v\\(v) = v—^-, (5.5) 
u E (v) = 2u a (v) - 2u D (v) - u\\(v), (5.6) 

where Erf(x) = e~ y2 dy is the error function, G(x) = [Erf(x)— xErf(x)]/2x 2 

is the Chandrasekhar function, and v = y/^nnQq^YtiAT^^m' 1 / 2 is the dimensional 
like-particle collision frequency (here In A is the Coulomb logarithm and q is the par- 
ticle charge). 

If one wishes to construct a model linearized collision operator, the following 
general form constitutes a natural starting point 



D[v) 'dvY 



+ P[Sf](v)F , (5.7) 



where the first term is the "test-particle" collision operator and the second term 
the "field-particle" operator. Most model operators can be obtained by picking a 
suitably simple form for the velocity-space diffusion tensor D and the functional P, 
subject to the constraints that one chooses to impose on the model operator. 

In constructing our model operator, we retain the exact form of D for the 
linearized Landau collision operator |67j : 



C[5f] = u D L[5f] + i|- QA^IQ + P[Sf](v)F , (5.8) 

where we have explicitly separated the energy-diffusion part (the second term) and 
the angular part (the first term), which includes pitch-angle scattering and is de- 
scribed by the Lorentz operator: 



L[Sf] - I 



d M , 2 ,d5f 1 d 2 5f 



dC d£ 1 -£ 2 d$ 2 

Our modeling choice is to pick P to be of the form 



(5.9) 



2v ■ U\Sf] v 2 
P[5f](v) = v s U 2 [ ° Jl + u E — Q[5f]. (5.10) 
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One can view this prescription as first expanding P in spherical harmonics (one can 
easily show that they are eigenf unctions of the full field-particle operator), retaining 
only the first two terms, and then arbitrarily factorizing the explicit v and Sf de- 
pendence of each harmonic. The functionals U[Sf] and Q[Sf] are mandated to have 
no explicit velocity dependence. In this ansatz the v dependence is chosen so that 
the final operator is self adjoint and also to ensure automatic particle conservation 
by the field-particle operator: J P[Sf](v)F dv = 0. Indeed the first term in Eq. 



(5.10) gives a vanishing contribution to this integral because it is proportional to v, 
and so does the second term because v^veFo = —(d/dv)(v 5 b , \\F ). The functionals 
U[5f] and Q[Sf] are now uniquely chosen so as to ensure that the model operator 
conserves momentum and energy: a straightforward calculation gives 



3 / v s vSf dv 
2 / (v/v th ) 2 u s F dv 



J v 2 (v/vth) vsFodv 

These are in fact just the standard correction expressions used for the model pitch- 
angle-scattering operator [681 EZ] and for more complex operators including energy 
diffusion [23] . 

To summarize, we now have the following model operator for like particle 
collisions: 



C[Sf] = v -f 



d 2 dSf 1 dHf 
dC 1 -i 2 dti 2 

+ Us ^wi Fo + UE ^ Q[6f]Fo , 



i d (i 4 ^ d Sf 



v 2 dv V 2 " dv F 



(5.13) 



where the functionals U[8f] and are given by Eqs. (5.11) and (5.12). The 

modeling choice of the field-particle operator that we have made [Eq. (5.10)] means 
that, in order to compute our collision operator, we have only to calculate definite 
integrals over the entirety of the velocity space — a significant simplification in 
terms of computational complexity and ease of use in numerical simulations (see 
Chapter 6). 

As we have shown above, our operator conserves particles, momentum and 
energy by construction. It is also not hard to see that it vanishes precisely when 
Sf/F = (1, v, v 2 ) and linear combinations thereof, i.e. if Sf is a perturbed Maxwellian. 
From this and the fact that the operator is self adjoint, it can be shown that the 
operator only conserves particles, momentum and energy and that no spurious con- 
servation laws have been introduced by our model. Because the operator contains 
the exact test-particle part, it provides velocity-space diffusion both in energy and 
in angle and thus will efficiently dissipate small-scale structure. Finally, it satisfies 
the if -theorem, as proved in Appendix [Cj 



Our operator thus fulfills the criteria set forth in Sec. |5.1| to be satisfied by 
a physically reasonable model operator. We now proceed to convert this operator 
into a form suitable for use in gyrokinetics. 
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5.3 Collisions in Gyrokinetics 



The gyrokinetic theory is traditionally derived for a collisionless plasma 0, EH] ■ 



However, as we have argued in Sec. 5.1 even when the collision frequency is small, 
collisions should be included in order to regularize the phase space and to ensure 
convergence of fluxes to statistically stationary values. Mathematically, collisions 
can be included in gyrokinetics if the collision frequency is formally ordered to be 
comparable to the fluctuation frequency [TO], v ~ uo ~ k^v^ — the weakly collisional 
limit (collisionality larger than this leads simply to fluid equations). In practice, 
collision frequency tends to be smaller than the fluctuation frequency, but this need 
not upset the formal ordering as long as it is not too small: the cases v ^> to and 
v <C oj can be treated as subsidiary limits [13] . 

Under the formal ordering v ~ oj (and, in fact, also under an even less restric- 
tive ordering allowing for even smaller collisions^]), we have shown in Chapter 2 that 
the equilibrium distribution function (lowest order in the gyrokinetic expansion) is 
a Maxwellian, and the full distribution function can be represented as 

/= (l-^j F + h(t,R,fi,e), (5.14) 

where F is a Maxwellian, <3> the electrostatic potential (a fluctuating quantity) and h 
the (perturbed) distribution function of the particle guiding centers. Here e = mv 2 /2 
is the particle energy, \i = mv\j2B^ the first adiabatic invariant, Bq the strength of 
the equilibrium magnetic field, -R = r — p = r — bxv/Q the guiding center position, 
Q the cyclotron frequency, and b = B /B . The gyrokinetic equation, written in 
general geometry and including the collision operator is then 

dh , - . dh c , . . , , 

= -i^§r a -^T + J- W, <x) H } + C GK [ft], (5.15) 







where x = & ~ v ■ A/c the gyrokinetic potential, (x)r = (V^ 71 ") / x{R + p) d& is 
an average over gyroangles holding R fixed (the "gyroaverage" ) , is the guiding 



center drift velocity defined in Eq. (3.41) of Chapter |3j 

The gyrokinetic collision operator Cqk[^] is the gyroaverage of the linearized 
collision operator. The latter acts on the perturbed distribution h holding the 
particle position r (not the guiding center R) fixed. The latter nuance must be 
kept in mind when working out the explicit form of Cgk[M-R)] from the unaveraged 
linearized operator C[h{r — p)]. 

Let us restrict our consideration to local simulations, which are carried out 
in a flux tube of long parallel extent, but short perpendicular extent. In such 
simulations, one assumes that the equilibrium profiles are constant across the tube, 
but have non-zero gradients across the tube so as to keep all the appropriate drifts 



2 S. C. Cowley, unpublished 
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and instabilities. This permits one to use periodic boundary conditions and perform 
the simulations spectrally perpendicular to field lines [IT] . Thus 



(5.16) 



where I is a coordinate along the field line and the Fourier transform is understood to 
be only with respect to the perpendicular components of R, i.e., k = k±. Treating 
the perpendicular coordinates spectrally confines all dependence on the gyroangle 
$ to the exponent, thus we can compute the gyroangle dependence explicitly and 
carry out the gyroaveraging of the collision operator in a particularly transparent 
analytical way [25j [15] : 



CgkM 




ikRi 



k 



(e ik ' r C[e- ik - p h k } 



R 



ik - R {e ik - p C[e- ikp h k \) 



R- 



R 



(5.17) 



where p = b x v±/Q. Thus, in Fourier space 

C GK [h h ] = {e ik - p C[e- ih -»h k )) , 



(5.18) 



where (...) refers to the explicit averaging over the $ dependence. Some general 
properties of this operator are discussed in Appendix B of Ref. 



We now apply the general gyroaveraging formula Eq. (5.18) to our model op 



erator given by Eq. (5.13). The gyrokinetic transformation of variables (r,v,£, $) — > 
(R, /i, e, i?) mixes position and velocity space. However, in the collision operator, 
to the lowest order in the gyrokinetic expansion, we can neglect spatial dependence 
of \x that comes via the equilibrium magnetic field B (r) and thus use the (v,£) 
velocity variables. After some straightforward algebra, which involves converting 
velocity derivatives at constant r to those at constant R and evaluating the arising 



gyroaverages as detailed in Appendix [Dj we arrive at the following model gyrokinetic 
collision operator 



2 di 



2 dhk + ]^d_ 
<9£ v 2 dv 



v 4 v\\F 



d_hk 

dvYn 



4 

+ 2v, 



[u D (i+e)+Mi-e 



k ± p h k 

v±Ji(a)U ± [hk\+v\\J (a)U\i[h k ] „ 



(5.19) 



'th 



F Q + v E -h- Jo{o)Q[h k \F Q 



y th "th 
where p = v^/Q is the thermal Larmor radius (not to be confused with the velocity- 
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dependent p), a = k±v±/Q, Jo and J\ are Bessel functions and 

U±[h k ] = ^ v s v±Ji(a)h k dv j J (v/v th ) 2 v s F dv, (5.20) 

U \\l h k\ = ^ J v s v\\J (a)h k dv j J (v/v th ) 2 u s F dv, (5.21) 

Q[hk) = v 2 u E J {a)h k dv / / v 2 {v/v t hf veFq dv. (5.22) 



Note that since the position and velocity space are mixed by the gyrokinetic trans- 
formation of variables, R = r — p, the collision operator now contains not just 
pitch-angle and v derivatives but also a spatial perpendicular "gyrodiffusion" term. 

It is important to make sure that the operator we have derived behaves in a 
physically sensible ways in the long- and short- wavelength limits. When k±p <C 1, 
all finite-Larmor-radius effects disappear, and we end up with pitch-angle scattering 
and energy diffusion corrected for energy and parallel momentum conservation - 
the drift-kinetic limit. In the opposite limit, k±p ^> 1, we can estimate the behavior 
of our operator by adopting the scaling of the velocity derivatives based on the non- 
linear perpendicular phase mixing mechanism for gyrokinetic turbulence proposed 
in Ref. [13] : this produces velocity-space structure with characteristic gradients 
Vthd/dv± ~ k±p (see also Refs. [60j [15]). With this estimate, we see that all the 
field-particle terms in the operator are subdominant by a factor of {k±p)~ 3 . Thus 
the operator reduces to the gyrokinetic form of the test-particle Landau operator 
in this limit. All diffusive terms are also equally large in this scaling, supporting 
our supposition that energy diffusion needs to be included. These considerations 
give us some confidence that we correctly model the diffusive aspects of the colli- 
sional physics in a short-wavelength turbulent regime. Indeed, if one applies the 
same estimates to the full linearized Landau operator, the Rosenbluth potentials of 
the perturbation are small when k±p ^> 1 because they are integrals of a rapidly 
oscillating function, so the dominant effect does, indeed, come entirely from the 
test-particle part of the operator. 



The gyrokinetic collision operator given by Eq. (5.19) respects the if-theorem: 

E 

jj Y o C GK [h]dRdv<0, (5.23) 

which is the what has to be satisfied in order for heating and transport in gy- 
rokinetics to be correctly calculated [TDl H5] . The operator also manifestly diffuses 
small-scale structure both in velocity and in (perpendicular) position space. One 
cannot, however, perform the conservation- law tests upon this operator because one 
cannot separate the position- and velocity-space dynamics, and hence collisional and 
collisionless dynamics, in the gyrokinetic phase space. Thus, we take the view that 
the conservation laws are guaranteed for the gyrokinetic collision operator in the 
sense that they were guaranteed for the original model operator from which it was 



This can either be shown directly from Eq. (5.191 (analogously to the proof in Appendix [c| 



or inferred from Eq. (5.2 I by transforming to gyrokinetic variables 
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derived. For practical numerical applications, this leaves the question of how this 
operator is best discretized and implemented. This is addressed in Chapter 6, where 
we also demonstrate the correct performance of our model operator on a number of 
test problems and show that all its new components (energy diffusion, gyrodiffusion, 
conservation terms) are necessary to avoid unphysical results. 



5.4 Electron-Ion Collisions 

Let us now turn to the collisions between different species and focus on a 
plasma containing only electrons and one species of ions with a mass ratio m e /mi <C 
1. The smallness of the mass ratio allows for a significant simplification of the 
interspecies collision terms. Since ion-electron collisions are subdominant to ion-ion 
collisions [HZ], v ie jvn ~ {m e jm^) 1 l 2 , it is safe for most physical purposes to neglect 
the ion-electron collisions and effects associated with them (such as the small slow 
collisional change in the mean ion momentum). Thus, the ion collisions can be 



modeled using the like-particle operator proposed above [Eq. (5.19) . 

The situation is different for the electron-ion collisions, which are the same 
order in mass ratio as the electron-electron collisions [HZ], v e i ~ v ee- Thus, the full 
electron collision operator has two parts: 

C[5f e ] = C ee [6f e ] + C ei [5f e \. (5.24) 

The electron-electron operator C ee [5f e ] can be modeled by the like-particle operator 



proposed above [Eq. (5.13)], the electrom-ion collision operator can be expanded in 



the mass ratio and to lowest order reads B 



2l? • u ■ 

C ei [5f e ] = u^[L[Sf e ] + -^F 0e ), (5.25) 



Vthe^ 3 



V the 

v e n(v) = Vei{^) (5.26) 

where z/ ei = \/27rn 0i Z 2 e 4: \n AT^^rrie 1 ^ 2 is the dimensional electron-ion collision 
frequency, Z = qi/e, e is the fundamental charge, L is the Lorentz operator given 
by Eq. (5.9), and 

1 



u t = — vSfrdv (5.27) 
n 0i J 

is the ion flow velocity. Thus, the electron-ion collisions are correctly modeled to 
lowest order in the mass ratio by electron pitch-angle scattering off static ions plus 
electron drag against the bulk ion flow. Note that the ion drag term is necessary 
to correctly capture electron-ion friction and hence resistivity; failure to include 
it leads to incorrect results, with mean electron momentum relaxed towards zero 
rather than towards equality with the mean ion momentum. 

Performing the conversion of C e i to the gyroaveraged form in a way analogous 



79 



to what was done in Sec. |5.3| and Appendix [DJ we get 



D 



2 dC 



dh 



ek 



V 



2 2 2v ll J (a e )u llik 

K±P e llek H o t{ 



Oe 



the 



^th e 



Zm e U| Ji(o e ) 2 2 1 
^Oefc_LPj 



mi U. 



the 



«0i 



2u', 2 JifoJ 



7 th« 



where 



1 

™0i 



ui|J (aj)/ij fe cfo 



(5.28) 



(5.29) 



and a s = k±v±/Q s for species s and the rest of the notation as the same as in 
previous sections, with species indices this time. 



Let us estimate the size of the four terms in Eq. (5.28) at the ion (long) and 



electron (short) scales. The first term (pitch-angle scattering) is always important. 
At the ion scales, kj_pi ~ 1, the third term (parallel ion drag) is equally important, 
while the second term (electron gyrodiffusion) and the fourth term are subdominant 
by a factor of m e /mj. At the electron scales, k±p e ~ 1, the pitch-angle scattering 
and the electron gyrodiffusion (the first two terms) are both important. Since at 
these scales k±pi ~ (mj/m e ) 1//2 ^> 1, the third and fourth terms are subdominant 
by a factor (resulting from the Bessel functions under the velocity integrals) of 
1 / \fk±p~i ~ (m e / nii) 1 ^. In fact, they are smaller than this estimate because at these 
short wavelengths, the ion distribution function has small-scale structure in velocity 
space with characteristics scales 5v±/vthi ~ l/k±pi, with leads to the reduction of 
the velocity integrals by another facto r of l/\Zk_ s _p i . Thus, at the electron scales, 
the third and fourth terms in Eq. (5.28) are subdominant by a factor of (m e /mi) l l 2 . 



These considerations mean that the fourth term in Eq. (5.28) is always negligi- 



ble and can safely be dropped. The full model gyrokinetic electron collision operator 
is therefore 



C GK [h 



efc 



D 



2x dh e k 



2 8C q J d£ 



1 



k\p\h ek + 



J th t 



2v» J (a e )u\\ ik 



Oe 



7 the 



(5.30) 



where the electron-electron model operator CQ K [h e k] is given by Eq. (5.19) and wiufe 
by Eq. KM. 



Finally, we note that since en 0e (u\\i — u\\ e ) = j\\ is the parallel current, the 
parallel Ampere's law can be used to express u\u in Eq. (5.30) in a form that does 



not contain an explicit dependence on the ion distribution function: 

If c 

/ Dii J (a e )h ek dv + 

n e J 4:7ien 0e 



U\\ik 



k±Au k . 



(5.31) 



This turns out to be useful in the numerical implementation of the electron operator, 
detailed in Chapter 6. 
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5.5 Summary 



In Sec. |5.1| we have argued the necessity of dissipation in turbulence simula- 
tions, justified the direct modeling of collisions in order to provide such dissipation 
and postulated a set of constraints for a physically reasonable model collision oper- 
ator. Previously used model operators were deemed unsatisfactory, in part because 
the majority of them do not contain a mechanism for energy diffusion. Two of 
the well-known existing model operators that contain energy diffusion are detailed 
in Refs. [25] and [21] . However, the former does not satisfy the if -theorem [Eq. 
(5.2)], and the latter incorrectly captures the smallest scales. These problems are 
demonstrated and discussed in detail in Appendix [Ej 



In Sec. 5.2 we presented a new operator [Eq. ( |5.13 )] that successfully intro- 
duces energy diffusion while maintaining the if-Theorem and conservation laws, 
thus satisfying the conditions set forth in the introduction. This operator is then 
transformed into gyrokinetic form in Sec. 5^3 correctly accounting for the gyrod- 
iffusive terms and FLR effects [Eq. ( |5.19 )]. In order to provide a complete recipe 
for modeling the collisional effects in simulations, the same gyroaveraging procedure 
is applied in Sec. |5.4| to electron-ion collisions, somewhat simplified by the mass- 
ratio expansion [Eq. (5.30)]. This leaves us with a complete picture of collisions 
in gyrokinetic simulations, capturing gyrodiffusion, resitivity and small-scale energy 
diffusion. 

When we discussed the gyroavergaing procedure in Sec. |5.3| we presented the 
specific case of the application to Eulerian flux-tube Sf gyrokinetic simulations [TD1 
153] . However, the form presented in Eq. (5.13) is suitable for inclusion in most 
Sf kinetic systems and even amenable to use in Lagrangian codes by applying the 
methods of Refs. [71] or [72J to the gyroaveraged operator given by Eq. (5.19). 
Indeed, by suitable discretization of the gyroaveraging procedure [55] it would also 
be usable in a global Eulerian code. 

We conclude by noting that the final arbiter of the practicality and effectiveness 
of this collision model is the numerical implementation and testing performed in 
Chapter 6, where our operator is integrated into the Trinity code. The battery 
of tests shows that our operator not only reproduces the correct physics in the 
weakly collisional regime but even allows a gyrokinetic code to capture correctly the 
collisional (reduced-MHD) limit. 
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Chapter 6 

Linearized model Fokker-Planck collision operator for 
gyrokinetics: numerics 



6.1 Introduction 

Collisions play an important role in gyrokinetics. An accurate collision opera- 
tor is important for calculation of neoclassical transport [73l [71] and the growth rate 
of instabilities such as trapped electron modes [T5J [76], dissipative drift waves [TTJ 
[THEE], and microtearing modes [79] in moderate collisionality regimes. Collisions 
can also affect the damping of zonal flows [SO] and other modes that provide a sink 
for turbulent energy. In their absence, arbitrarily fine scales can develop in phase 
space jHU [53j [15], El El] , which can in some cases pose challenges for discrete nu- 
merical algorithms, especially in the long-time limit [821 183] : even a modest amount 
of collisions can make accurate numerical calculation much easier. 

Furthermore, inclusion of a small collisionality keeps the distribution function 
smooth enough in velocity space that the standard gyrokinetic ordering [9] for ve- 
locity space gradients is satisfied. For example, the parallel nonlinearity [8^ 185] . 
given by 

_ JL 

dv { \ 

enters at the same order as the other terms in the gyrokinetic equation if the typical 
scale of parallel velocity fluctuations, 5v\\, is one order smaller in the gyrokinetic 
expansion parameter p/L (p = gyroradius and L = background scale length) than 
the thermal speed, v t h- Here, h is the non-Boltzmann part of the perturbed distri- 
bution function (defined more rigorously in the next section), $ is the electrostatic 
potential, B is the magnetic field strength, b = H /B , q is particle charge, m is 
particle mass, and ( . ) denotes the gyroaverage at fixed guiding center position R. 

While such a situation is possible in the collisionless limit, a small collisionality 
prohibits the formation of structures with Sv» ~ {p/L)v t h- The level of collisionality 
necessary to negate the importance of the parallel nonlinearity can be calculated by 
assuming a balance between collisions and fluctation dynamics: 

f - C[h] =► a* ~ (6.2) 
where C[h] describes the effect of collisions on h, uo is the fluctuation frequency, 



, . q - b x VB . ^ , T , 



(6.1) 
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and v is the collision frequency. From the above expression, we see that scales in 
velocity space become small enough for the parallel nonlinearity to be important 
only when the collision frequency satisfies v ~ {p/L) 2 uo. Such low collisionalities 
are not present in most fusion plasmas of interest. Furthermore, if such an ordering 
had to be adopted, the lowest order distribution function could become strongly 
non-Maxwellian. This is clearly a problem for 5f codes that assume an equilibrium 
Maxwellian. 

In light of the above considerations, it is important to include an accurate 
treatment of dissipation in gyrokinetic simulations. In order to faithfully represent 
gyrokinetic plasma dynamics at reasonable numerical expense, we take the view that 
the form of the dissipation should be such that it: ensures satisfaction of the standard 
gyrokinetic ordering; locally conserves particle number, momentum, and energy; 
satisfies Boltzmann's if-Theorem; and efficiently smooths phase space structure. 
The first of these requirements has already been discussed in the context of the 
parallel nonlinearity. Conservation properties have been found to be important, for 
instance, in calculations of the neoclassical ion thermal conductivity (29], as well as 
in a wide range of problems in fluid dynamics. The existence of an if-Theorem is 
critical for entropy balance [5^ 1621 ITU] and for the dynamics of the turbulent phase 
space cascade [T51 [J3]. Efficient smoothing of phase space structures is necessary to 
resolve numerical simulations at reasonable computational expense. 

A commonly employed dissipation mechanism in gyrokinetic simulations is ar- 
tificial (hyper) dissipation, often in physical (position) space [51j [7UJ |5UJ ED 186] . 
Ideally, the form of the artificial dissipation should be chosen to satisfy the require- 
ments listed above and should be tested for convergence to the collisionless result. 
Of course, artificial dissipation alone is unable to capture the correct dynamics for 
moderate to strongly collisional systems where turbulent fluxes and other observable 
quantities depend sensitively on collisionality; for such systems, a physical dissipa- 
tion model is desired. 

A number of such model physical collision operators are employed in gyroki- 
netic codes [m [701 EH EZ] ■ These range in complexity from the Krook operator [88] 
to the Rutherford-Kovrizhnikh operator [68J to the Catto-Tsang operator [25J, all 
of which have previously been implemented in GS2 (see, e.g. Ref. |87j). However, 
none of these satisfy all of the properties we require of a good collision operator 
(See Appendix |E] for a fuller discussion of this point). Here, we discuss numerical 
implementation in Trinity of an improved model operator which: includes the ef- 
fects of both pitch-angle scattering and energy diffusion (i.e. efficiently smooths in 
phase space and ensures gyrokinetic ordering); conserves particle number, momen- 
tum, and energy; satisfies Boltzmann's if-Theorem; and reduces to the linearized 
Landau test-particle operator in the large k±p limit. A full description of this op- 
erator and a discussion of its desirable properties is given in Chapter 5. We will 
focus on how such an operator can be implemented efficiently in gyrokinetic codes 
while maintaining the properties listed above and on how our gyrokinetic dissipation 
scheme (or any other) might be tested against a number of plasma physics problems. 
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This chapteiQis organized as follows: in Sec. (12 we present the gyroaveraged 
collision operator derived in Chapter 5 and examine properties that should be taken 
into account when using it in numerical simulations; in Sec. 6.3, we describe our 



numerical implementation of the collision operator; in Sec. 6.4 we present numerical 
results for a number of tests demonstrating the ability of our collision operator 
implementation to reproduce correct collisional and collisionless physics; and in Sec. 
6.5[ we summarize our findings. 



6.2 Properties of the gyroaveraged collision 
operator 



In order to include collisions in gyrokinetics, we follow the treatment of Ref. [10] 
and assume the collision frequency, u, to be the same order in the gyrokinetic or- 
dering as the characteristic fluctuation frequency, As was shown in Chapter 2, 
this leads to the requirement that the distribution of particles in velocity space is 
Maxwellian to lowest order and allows us to represent the total distribution function 
through first order in p/L (where p is ion gyroradius and L is the scale length of 
equilibrium quantities) as 

/(r, p, e, t) = F (e) (l - ^|^J + h(R, p, e, t), (6.3) 

where r is particle position, R = r — b x v/Qq is guiding center position, p = 
mv'j_/2Bo is magnetic moment, e = mv 2 /2 is particle energy, Fq is a Maxwellian, $ 
is the electrostatic potential, B is the magnitude of the background magnetic field, 
Tq is the background temperature, q is particle charge, and fl = qB /mc. The 
gyrokinetic equation governing the evolution of h is given by 

Oh / \ dh c 

S + ("i b + VD )-aa + Bb {W - k} (64) 

where b = Hq/Bq, vb is the drift velocity of guiding centers, x = $ ~ v ' A/c, A is 
the vector potential, {a, b} is the Poisson bracket of a and b, (a) R is the gyroaverage 
of a at constant R, and (C[/i]) R is the gyroaveraged collision operator. 

For (C[h]) n , we restrict our attention to the model collision operator presented 
in Chapter 5. We work within the framework of the continuum gyrokinetic code 
GS2 [51J, which assumes periodicity in the spatial directions perpendicular to B 
in order to reduce the simulation volume to a thin flux tube encompassing a single 
magnetic field line. Consequently, we require a spectral representation of (C[h]) R : 

(C[h)) n = J2 eik ' nC G«l h ^ ( 6 ' 5 ) 
k 

1 This chapter taken from Ref. [22] , 

2 Note that this ordering does not prevent one from considering the cases of v -c u> and v ^S> lo 
as subsidiary orderings [T3] 
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where k is the perpendicular wavevector. For convenience, we reproduce the expres- 
sion for the same-species part of CcK[h k ] from Chapter 5 in operator form: 



where 



and 



C GK [K] = L[h k ] + D[h k ] + U L [h k ] + U D [h k ] + E[h k ], (6.6) 

L[h k ) ^^(l-e) 9 ^- ^vn (1 + e) K (6.7) 
r n I d f A d h k \ k 2 v 2 / ox , / % 

D ^ s 2^ h« ^ ) - «r» f 1 - « > ftk (6 - 8) 

are the gyroaveraged Lorentz and energy diffusion operators (which together form 
the test-particle piece of the linearized Landau operator, as shown in Refs. [25] 
and [21]), 

tt \h i - t? ( t ( \ I d?vv D v\\Jo{a)K J (Pv iy D v ± J 1 (a)h k \ 
U L [h k \ = u D F [J {a)v\\ tt— + Ja(a)v ± jr- ^- j 6.9 



U D [h k \ = -AuF [ J a f || , — ^ + J i (°) U -L r « a — 2^ ) 



and 

J G? 3 t> Ai/uy Jo(a)h k J d 3 v Auv±Ji(a)h k s 

i 

(6.10) 

are the gyroaveraged momentum-conserving corrections to the Lorentz and energy 
diffusion operators, and 

E[h k \ = v E v J (a)F tt— — 6.11 

J drv v E i ' F 

is the gyroaveraged energy-conserving correction (the conserving terms are an ap- 
proximation to the field-particle piece of the linearized Landau operator). The elec- 
tron collision operator has the following additional term to account for electron-ion 
collisions: 



\h i - u ei ( l 9 (i P\ dheM kW h +r 2 \h 4. 2v w u w^ u„\w\ 

^GKi^eM - p d[^-q? I 1 -? ) - + £ ) ^,k+ — 2 J {a e )£ 0e j, 

^ Oe th e 

(6.12) 

where £ = v\\/v is the pitch angle, a = kv±/Qo, Jo and J\ are Bessel functions of 
the first kind, v t h = a/ 2T / m is the thermal velocity, and u\\ [h i:k ] is the perturbed 
parallel ion flow velocity. Expressions for the velocity-dependent collision frequencies 
up, Au, un, and Ve are given in Chapter 5 (which follows the notation of Ref. [21]). 

Having specified the form of our collision operator, we now discuss some of its 
fundamental properties that guide our choice of numerical implementation. 
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6.2.1 Collision operator amplitude 



Even when the collisionality approaches zero, C G k\)i]s\ can have appreciable 
amplitude. There are two reasons for this: first, the velocity dependence of u D , u E , 
and Az/ is such that each go to infinity as v — > (so low-velocity particles are always 
collisional); and second, we expect the distribution function to develop increasingly 
smaller scales in v and £ as collisionality decreases, so that the amplitude of the 
terms proportional to d 2 h/d£ 2 and d 2 h/dv 2 may remain approximately constant]^] 
(i.e. C GK [h k ] y4 as v -> 0) [Ml G2H E51 EQ. The fact that C GK [h k ) can be quite 
large even at very low collisionalities means that it should be treated implicitly if 
one wants to avoid a stability limit on the size of the time step, At. In Sec. |6.3[ we 
describe our fully implicit implementation of the collision operator. 



6.2.2 Local moment conservation 

Since collisions locally conserve particle density, momentum, and energy, one 
would like these properties to be guaranteed by the discrete version of the collision 
operator. Mathematically, this means that the density, momentum, and energy 
moments of the original (un-gyroaveraged) collision operator must vanish (for same- 
species collisions). However, the non-local nature of the gyroaveraging operation 
introduces finite Larmor radius (FLR) effects that lead to nonzero values for the 
analogous moments of (C[/i]) R . Since this is the quantity we employ in gyrokinetics, 
we need to find the pertinent relations its moments must satisfy in order to guarantee 
local conservation properties. 

This is accomplished by Taylor expanding the Bessel functions Jo and J\. In 
particular, one can show that [21] 



/ 




d% | v | ((C[h]) R ) F = ^ e* kr / d*v [ Vll b\ C GK [K) - V • T C , (6.13) 



where (.) r denotes a gyroaverage at fixed r, C GK [h]i\ is the operator of Eq. (6.6 ) with 
kp = (neglecting FLR terms but retaining nonzero subscripts k for h), and Tc is 
the collisional flux of number, momentum, and energy arising from FLR terms. Con- 
sequently, the density, momentum, and energy moments of the gyrokinetic equation 
can be written in the conservative form 

d ^ + V-Y M = I d 3 v[ Vll b\ C GK [K], (6.14) 

where Ai = (5n nSu\\ 8p) T represents the perturbed number, momentum, and en- 
ergy densities, the superscript T denotes the transpose, and Tm contains both the 
collisional flux, Tq, and the flux arising from all other terms in the gyrokinetic 




3 This is analagous to the result in fluid turbulence where the dissipation rate remains finite as 
viscosity becomes vanishingly small. 
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equation (for a more detailed discussion, see Ref. [21] ). Thus, local conservation 
properties are assured in gyrokinetics as long as the density, momentum, and en- 
ergy moments of C^^k] vanish: 

J d 3 v ^jj C° GK [h k ] = 0. (6.15) 



We describe how this is accomplished numerically in Sec. 6.3 



6.2.3 //-Theorem 

In contrast with local conservation properties, the statement of the if -Theorem 
is unmodified by gyroaveraging the collision operator. Defining the entropy as S = 
—/In/, Boltzmann's if -Theorem tells us 



as rd 3 



j^jd'v ln[/]C[/]>0, (6.16) 



where V = J d 3 r and the double integration spans phase space (the velocity integra- 
tion is taken at constant particle position r). Expanding the distribution function 
as before, we find to lowest order in the gyrokinetic ordering 

^Jd 3 v^C[h}<0. (6.17) 
Changing variables from particle position r to guiding center position R, we obtain 

d3v J^Y <CW>R < 0, (6.18) 

where now the velocity integration is taken at constant R. In this case, the non- 
locality of the gyroaveraging operation leads to no modification of the ii-Theorem 
because of the definition of entropy as a phase-space averaged quantity (as opposed 
to local conservation properties, which involve only velocity-space averages). There- 
fore, one can easily diagnose entropy generation and test numerical satisfaction of 



the ii-Theorem in gyrokinetic simulations, as we show in Sec. 6.4 



6.3 Numerical implementation 

It is convenient for numerical purposes to separately treat collisional and col- 
lisionless physics. Thus, we begin by writing the gyrokinetic equation in the form 

r)h 

-^ = C GK [h k ]+A[h k }, (6.19) 
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where A[hk] represents the rate of change of due to the collisionless physics. In 
order to separate these terms, we utilize Godunov dimensional splitting [89], which 
is accurate to first order in the timestep At: 

h ^^ = AK,K\ (6.20) 
h n+i — h* 

*±^ = C aK [hp\ (6.21) 

where n and n + 1 are indices representing the current and future time steps, and 
/t£ is defined by Eq. (6.20) - it is the result of advancing the collisionless part of 
the gyrokinetic equation. With /i£ thus given, we restrict our attention to solving 
Eq. (6.21). For notational convenience, we suppress all further k subscripts, as we 
will be working exclusively in fc-space. 

As argued in Sec. 6.2 , we must treat the collision operator implicitly to avoid 
a stability limit on the size of At. We use a first order accurate backward- difference 
scheme in time instead of a second order scheme (such as Crank-Nicholson |90j) 
because it is well known that the Crank-Nicholson scheme introduces spurious be- 
havior in solutions to diffusion equations when taking large timesteps (and because 
Godunov splitting is only first order accurate for multiple splittings, which will be 
introduced shortly). 

With this choice, h n+l is given by 

h n+1 = (1 - AtCcK f 1 h*. (6.22) 

In general, Cqk is a dense matrix, with both energy and pitch-angle indices. Inver- 
sion of such a matrix, which is necessary to solve for h n+l in our implicit scheme, 
is computationally expensive. We avoid this by taking two additional simplifying 
steps. First we employ another application of the Godunov splitting technique, 
which, combined with the choice of a (£, v) grid in Trinity [81], allows us to con- 
sider energy and pitch-angle dependence separately: ^ 

h** = [1 — At (L + U L )Y l h* (6.23) 
h n+1 = [1 - AtiD + Uo + E^h**, (6.24) 

The h* and h** are vectors whose components are the values of h at each of the 
(£, v) grid points. In Eq. ( |6.23[ ) we order the components so that 

h = (/in, h 2 \, hm, hu, ■-, h NM ) T , (6.25) 

where the first index represents pitch-angle, the second represents energy, and iV 
and M are the number of pitch-angle and energy grid points, respectively. This 
allows for a compact representation in pitch-angle. When solving Eq. (6.24), we 
reorder the components of h so that 

h = (h u , h 12 , hi N , h 2 i, h NM ) T , (6.26) 

allowing for a compact representation in energy. 



4 We note that an implicit treatment of the Catto-Tsang operator (including energy diffusion) 
has independently been implemented in GS2 using the same splitting technique [87 . 
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6.3.1 Conserving terms 



The matrices 1 — AtL and 1 — AtD are chosen to be tridiagonal by employing 
three-point stencils for finite differencing in £ and v. This permits computationally 
inexpensive matrix inversion. However, the full matrices to be inverted include the 
momentum- and energy-conserving operators, U and E, which are dense matrices. 
We avoid direct inversion of these matrices by employing the Sherman-Morrison 
formula E2], which gives x in the matrix equation Mx = b, as long as M can 
be written in the following form: 



M = A + u <g> v, 
where Cg) is the tensor product. The solution is then given by 

v-y 



V ■ z 



(6.27) 



(6.28) 



where y = A _1 b, z = A _1 u, and the dot products represent integrals over velocity 
space. If A -1 is known or easily obtainable (as in our case), this formulation pro- 
vides significant computational savings over the straighforward method of directly 
inverting the dense matrix M. 

Details of the application of the Sherman-Morrison formula to Eqs 



(6.23) 



and (6.24) are given in Appendix [FJ Here, we state the main points. The matrix 
operators L and D are to be identified with A, and the integral conserving terms U 
and E can be written in the form of the tensor product, u £g> v. Identifying h** and 
h n+1 with x, we find that multiple applications of the Sherman-Morrison formula 
give 

v 2 • y 2 



h 



Y2 



v 2 • z 2 



Z2, 



where 



Y2 = yo 

z 2 = z 



v • yo 


s - 


vi • y 


1 + v ■ s _ 


1 + Vi ■ w _ 


v • z 


s - 


Vl • z 


1 + v • s _ 


1 + Vi ■ w _ 



w , 
w . 



(6.29) 

(6.30) 
(6.31) 



The quantities v , Vl, v 2 , z , s , w , and y are specified in Table [Rl] in Appendix 
[Fj With the exception of y , each of these quantities is time- independent, so they 
need be computed only once at the beginning of each simulation. Consequently, 
inclusion of the conserving terms in our implicit scheme comes at little additional 
expense. 



We note that when Eq. (6.29) is applied to computing the inverse matrix in 



Eq. (6.23), the corresponding v 2 is nonzero only for the electron collision operator. 



This term arises by using the parallel component of Ampere's law to rewrite the 
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electron-ion collison operator of Eq. (6.12) as 



Cq K [h e 



1 d 



k 2 v 2 



(i + £ 2 ) K 



+ 



2vn 



-Jo(a e )F 0e 



J th e 



(6.32) 



where e is the magnitude of the electron charge, un [h e ] is the parallel component of 
the electron fluid velocity, and n e is the equilibrium electron density. For electron 
collisions, the A\\ term is absorbed into h* so that we use the modified quantity 



K = K + v*At 



ck 2 v\\ 



27rev4 n 0) , 



-A\\ J (a e )F 0jl 



(6.33) 



when applying the Sherman-Morrison formula, where An from the n + 1 time level 
is used (for details on the implicit calculation of Au, see Ref. [5T]). 



6.3.2 Discretization in energy and pitch angle 

We still must specify our choice of discretization for Cqk- Ideally, we would 
like the discrete scheme to guarantee the conservation properties and if -Theorem 
associated with C. As discussed in Sec. 6^2, the former is equivalent to requiring 
that the kp = component of Cqk, Cq K , satisfy Eq. (6.15). We now proceed to 
show that this requirement is satisfied by carefully discretizing the conserving terms 
and by employing a novel finite difference scheme that incorporates the weights 
associated with our numerical integration scheme. 

We begin by writing Cq K [K\ for same-species collisions: 



C GK [h] 



vo_d_ 
Aw \\F C 



dh 



J d 3 v Avvuh 
J d 3 v AuvjF c 



J__d_ 
2v 2 dv 



d h 



V\\V ^0 



+ u E v 2 F 



dvF 
J d 3 v v E v 2 h 



+ V D V\\Fq 



f d 3 v u D v\\h 
J d 3 v iy D VnF 



Id 3 



V V E V A Fq 



(6.34) 



With Cq K thus specified, we now consider numerical evaluation of the relevant 



moments of Eq. (6.34). To satisfy number conservation (J d 3 v Cq K [}i] = 0), velocity 



space integrals of each of the terms in Eq. (6.34) should vanish individually. For 



integrals of the first two terms to vanish, we require a finite difference scheme that 
satisfies a discrete analog of the Fundamental Theorem of Calculus (i.e. conservative 
differencing); for the last three terms, we must have a discrete integration scheme 
satisfying f d 3 v u D v\\F = f d 3 v Avv\\F = f d 3 v u E v 2 F = 0. The requirement 
that f d 3 v z/ D t>||F = f d 3 v Az/t>||F = is satisfied by any integration scheme 
with velocity space grid points and associated integration weights symmetric about 
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v\\ = 0, which is true for the (£, v) grid used in Trinity. By substituting for Ve 
everywhere using the identity 



2 17 

v E v F 



1_ d_ 

v 2 dv 



(6.35) 



the other integral constraint (J d 3 v u e v 2 Fq = 0) reduces to the requirement that 
finite difference schemes must satisfy the Fundamental Theorem of Calculus. 

Parallel momentum conservation (J d 3 v v»CQ K [h] = 0) introduces the addi- 
tional requirements that: 



and 



d 3 v 



d v v\\ 



v || d 
2v 2 dv 



2 d£ v ' o£ J d 6 v v D vf,Fa 







4 17 
V\\V F 



d h 
dvY 



d 3 v Auv 2 F 



J d 3 v Ai/v\\h 
'jd 3 v Auv 2 F ' 



(6.36) 



(6.37) 



If the finite difference scheme used for all differentiation possesses a discrete version 
of integration by parts (upon double application), then Eqs. (6.36) and (6.37) are 
numerically satisfied as long as: v\\u D h in the second term of Eq. (6.36) is expressed 
in the form 



V\\VY)h = — - 



0_ 

2 V.^ 



dv\\ 



u D h, 



Av on the righthand side of Eq. (6.37) is expressed using the identity 

dv I " dv 



2Avv 3 Fn 



(6.38) 



(6.39) 



and all integrals are computed using the same numerical integration scheme (if 
analytic results for the integral denominators in terms three and four of Eq. (6.34) 
are used, then the necessary exact cancellation in Eqs. (6.36) and (6.37) will not 
occur). 

The only additional constraint imposed by the energy conservation require- 
GK [h] = 0) is that the form of Eq. (6.35) be slightly modified so 



ment (J d 3 
that 



v v 2 C% 



v E v F 



v 2 dv 



v\\v 4 F - 



dv 2 
dv 



(6.40) 



which still satisfies the number conservation contraint. Using the forms given by 
Eqs. (6.38)-(6.40), conservation properties are guaranteed as long as one employs a 
finite difference scheme for pitch-angle scattering and energy diffusion that satisfies 
discrete versions of the Fundamental Theorem of Calculus and integration by parts. 

For the case of equally spaced grid points in v and £, there is a straightforward 
difference scheme, accurate to secord order in the grid spacing, that satisfies both 
requirements [93] : 



—G— 

dx dx 



Gj+1/2 {hj+i ~ hj) — Gj-1/2 (hj — hj-i) 



Ax 2 



(6.41) 
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Number of grid points Number of grid points 



Figure 6.1: (Left): Solid line indicates the scaling of the leading order error, averaged 
over all grid points, of the conservative finite difference scheme for a Gauss-Legendre 
grid (the grid used in Trinity). The slope of the dotted line corresponds to a first 
order scheme. (Right): factor by which the conservative finite difference scheme of 
Eq. (6.42) amplifies the true collision operator amplitude at the boundaries of the 
Gauss-Legendre grid. 



where x is a dummy variable representing either v or £, Ax is the grid spacing, hj 
is the value of h evaluated at the grid point Xj, Xj±y 2 = (xj + Xj±x)/2, and G is 
either 1 — £ 2 (for pitch-angle scattering) or v^v^F® (for energy diffusion). However, 
in order to achieve higer order accuracy in the calculation of the velocity space 
integrals necessary to obtain electromagnetic fields, GS2 [81 J and a number of other 
gyrokinetic codes [55] use grids with unequal spacing in v and £ and integration 
weights that are not equal to the grid spacings. 

Given the constraints of a three-point stencil on an unequally spaced grid, we 
are forced to choose between a higher order scheme (a second order accurate scheme 



can be obtained with compact differencing (94], as described in Appendix G) that 
does not satisfy our two requirements and a lower order scheme that does. Since 
our analytic expression for Gqk was designed in large part to satisfy conservation 
properties (and because the conserving terms are only a zeroth order accurate ap- 
proximation to the field-particle piece of the linearized Landau operator [24"]). we 
choose the lower order scheme, given here, as the default: 

—G— « — (g ! . z hj+1 ~ k] - G ■_! , (6.42) 

dx dx Wj \ ^ +l ^ 2 Xj + \ — x,j * Xj — Xj-i J ' 

where Wj is the integration weight associated with Xj. 

Defining ^ = Gh', with the prime denoting differentiation with respect to x, 
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Figure 6.2: Plots showing evolution of the perturbed local density, parallel momen- 
tum, and energy over fifty collision times. Without the conserving terms (6.9 )-(6.1 1 ), 
both parallel momentum and energy decay significantly over a few collision times 
(long dashed lines). Inclusion of conserving terms with the conservative scheme 
detailed in Sec. 6.3 leads to exact moment conservation (solid lines). Use of a 
non-conservative scheme leads to inexact conservation that depends on grid spacing 
(short dashed lines). 
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Taylor series can be used to show 



(6.43) 



where Axj = Xj+1/2 — Xj-i/2- With the exception of pitch angles corresponding to 
trapped particles [EH EI], the grid points {xj} and associated integration weights 
{wj} in Trinity are chosen according to Gauss-Legendre quadrature rules [56]. For 



this case, we show numerically in Fig. 6.1 that 



1 N 

-Y 



AXn 



W ; 



1 + 




1 + 



and 



max 

=2,...,iV-l 











V iv) 


Wj 







(6.44) 



(6.45) 



where AT is the number of grid points in x. 

The boundary points (j = 1,N) are excluded from the max operator above. 



This is because Ax/w (the factor multiplying \E^- in Eq. (6.43 )) converges to approxi- 
mately 1.2 for the boundary points as the grid spacing is decreased (Fig. 6.1[ ). For the 
energy diffusion operator, we can make use of the property that G(x) = G(x)' = 
at x = and x = 00 to show 



± 



j±l/2 



% + o 



'(Ax) 



(6.46) 



with the plus sign corresponding to j = 1 and the minus sign to j = N. This is not 
true for the Lorentz operator, so we are forced to accept an approximately twenty 
percent magnification of the Lorentz operator amplitude at £ = ±1 and at the 
trapped-passing boundaries. We find that this relatively small error at the bound- 
aries has a negligible effect on measureable (velocity space averaged) quantities in 
our simulations. 



6.4 Numerical tests 

We now proceed to demonstrate the validity of our collision operator imple- 
mentation. In particular, we demonstrate conservation properties, satisfaction of 
Boltzmann's if- Theorem, efficient smoothing in velocity space, and recovery of the- 
oretically expected results in both collisional (fluid) and collisionless limits. While 
we do not claim that the suite of tests we have performed is exhaustive, it constitutes 
a convenient set of numerical benchmarks that can be used for validating collision 
operators in gyrokinetics. 
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6.4.1 Homogeneous plasma slab 



We first consider the long wavelength limit of a homogeneous plasma slab with 
Boltzmann electrons and no variation along the background magnetic field (k\\ = 0). 
The gyrokinetic equation for this system simplifies to 

d(6f) 



dt 



C° GK [h], 



(6.47) 



which means local density, momentum, and energy should be conserved. In Fig. |6.2 



we show numerical results for the time evolution of the local density, momentum, 
and energy for this system. 

Without inclusion of the conserving terms ( 6.9 )-(6. 1 1 ), we see that density is 
conserved, as guaranteed by the conservative differencing scheme, while the momen- 
tum and energy decay away over several collision times. Inclusion of the conserving 
terms provides us with exact (up to numerical precision) conservation of number, 
momentum, and energy. To illustrate the utility of our conservative implementa- 
tion, we also present results from a numerical scheme that does not make use of 
Eqs. (6.38)-(6.40) and that employs a finite difference scheme that does not possess 
discrete versions of the Fundamental Theorem of Calculus and integration by parts. 
Specifically, we consider a first order accurate finite difference scheme similar to that 
given by Eq. (6.42), with the only difference being that the weights in the denom- 
inator are replaced with the local grid spacings. In this case, we see that density, 
momentum, and energy are not exactly conserved (how well they are conserved de- 
pends on velocity space resolution, which is 16 pitch angles and 16 energies for the 
run considered here). 

The rate at which our collision operator generates entropy in the homogenous 
plasma slab is shown in Fig. CL3 As required by the if-Theorem, the rate of entropy 
production is always nonnegative and approaches zero in the long-time limit as 
the distribution function approaches a shifted Maxwellian. We find this to hold 
independent of both the grid spacing in velocity space and the initial condition for 
the distribution function (in Fig. 6.3, the values of h(£,v) were drawn randomly 
from the uniform distribution on the interval [—1/2, 1/2]). 



6.4.2 Resistive damping 

We now modify the system above by adding a finite A» . From fluid theory, we 
know that collisional friction between electrons and ions provides resistivity which 
leads to the decay of current profiles. Because the resistive time is long compared 
to the collision time, one can neglect d(Sf)/dt. However, since An ~ k~ 2 , and we 
are considering < 1, dA\\/dt must be retained. The resulting electron equation is 
of the form of the classical Spitzer problem (see, e.g., Ref. [67]): 

The parallel current evolution for this system is given by 

J|,(£) = J\\{t = 0)e^\ (6.49) 
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Figure 6.3: Plot of the evolution of entropy generation for the homogeneous plasma 
slab over twenty collision times. Our initial distribution in velocity space is random 
noise, and we use a grid with 16 pitch angles and 8 energies. The entropy generation 
rate is always nonnegative and approaches zero in the long-time limit. 

where 77 = I/um is the resistivity, an = 1.98r e n e e 2 /m e c 2 is the Spitzer conductivity, 
and r e = 3y / 7r/4z/ ei is the electron collision time. 

We demonstrate that the numerical implementation of our operator correctly 
captures this resistive damping in Figs. |6.4 and 6.5 We also see in these figures that 



in the absence of the ion drag term from Eq. (6.12), the electron flow is incorrectly 



damped to zero (instead of to the ion flow), leading to a steady-state current. 
6.4.3 Slow mode damping 

We next consider the damping of the slow mode in a homogenous plasma 
slab as a function of collisionality. In the low k±pi, high $ limit, one can obtain 
analytic expressions for the damping rate in both the collisional (kn\ m f p <C 1) and 
collisionless (knX m f p ^> 1) regimes, where X m f P is the ion mean free path (see e.g. 
Ref. [15J). The expressions are 



for k\\\ m f p <C 1, and 



<«■«» 



to = — r ' - (6.51) 
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Figure 6.4: Evolution of |Jii| for the electromagnetic plasma slab with (3 = 10~ 4 , 
kypi = 0.1, and v ei = \Qk\\v t h,i- Inclusion of the ion drag term in the electron-ion 
collision operator leads to the theoretically predicted damping rate for the parallel 
current given in Eq. (6.49) . Without the ion drag term, the parallel current decays 
past zero (at t ~ 22) and converges to a negative value as the electron flow damps 
to zero. 



97 



io- 3 



ion \ 

electron 

ion, no drag \ 

electron, no drag 

_i i i i I i i i i I i i i i_ 



10 20 30 

Time [(k„v t .)-«] 



Figure 6.5: Evolution of perturbed parallel flow for the electromagnetic plasma slab 
with (3 = 10~ 4 , kyPi = 0.1, and v ei = lOk^vu. Without inclusion of the ion drag 
term in Eq. (6.12), the electron flow is erroneously damped to zero (instead of to 
the ion flow). 



for k\\X m f p 3> 1. Here, va = v t h,i/VWi is the Alfven speed, and z/y is the parallel 
ion viscosity, which is inversely proportial to the ion-ion collision frequency, va\ 
u \\,i oc Vth.i/ vn- As one would expect, the damping in the strongly collisional regime 
[Eq. (6.50)] is due primarily to viscosity, while the collisionless regime [Eq. (6.51) 
is dominated by Barnes damping 



In Fig. |6.6[ we plot the collisional dependence of the damping rate of the 
slow mode obtained numerically using the new collision operator implementation 
in Trinity. In order to isolate the slow mode in these simulations, we took $ = 
An = 5n e = and measured the damping rate of 5B». This is possible because 5B\\ 
effectively decouples from $ and Au for our system, and 5n e can be neglected because 
Pi 3> 1 |15j. We find quantitative agreement with the analytic expressions (6.50) 
and (6.51) in the appropriate regimes. In particular, we recover the correct viscous 
behavior in the k\\X m f p <C 1 limit (damping rate proportional to v\\ t i), the correct 
collisional damping in the k\\X m f p ~ 1 limit (damping rate inversely proportial to 
v\\,i), and the correct collisionless (i.e. Barnes) damping in the k\\X m f p ^> 1 limit. 



6.4.4 Electrostatic turbulence 

Finally, we illustrate the utility of our collision operator in a nonlinear simula- 
tion of electrostatic turbulence in a Z-pinch field configuration [10] . We consider the 
Z-pinch because it contains much of the physics of toroidal configurations (i.e. cur- 
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Figure 6.6: Damping rate of the slow mode for a range of collisionalities spanning 
the collisionless to strongly collisional regimes. Dashed lines correspond to the 
theoretical prediction for the damping rate in the collisional (k\\X m f p <C 1) and 
collisionless (ku\ m f p ^> 1) limits. The solid line is the result obtained numerically 
with Trinity. Vertical dotted lines denote approximate regions (collisional and 
collisionless) for which the analytic theory is valid. 
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Figure 6.7: Evolution of ion particle and heat fluxes for an electrostatic, 2-species 
Z-pinch simulation. We are considering R/L n = 2.0 and v„ = 0.01v t h,i/ R- The 
particle flux is indicated by the solid line and is given in units of (p/R)no t iV t h,i- The 
heat flux is indicated by the dashed line and is given in units of (p/R)n i vf hi . We 
see that a steady-state is achieved for both fluxes without artificial dissipation. 



vature) without some of the complexity (no particle trapping). At relatively weak 
pressure gradients, the dominant gyrokinetic linear instability in the Z-pinch is the 
entropy mode [321 EU ESI ESI EZ1 EH], which is nonlinearly unstable to secondary 
instabilities such as Kelvin- Helmholtz |96j . 

In previous numerical investigations of linear [38] and nonlinear [96] plasma 
dynamics in a Z-pinch, collisions were found to play an important role in the damping 
of zonal flows and in providing an effective energy cutoff at short wavelengths. 
However, as pointed out in Ref. [HS], the Lorentz collision operator used in those 
investigations provided insufficient damping of short wavelength structures to obtain 
steady-state fluxes. Consequently, a model hyper- viscosity had to be employed. 

We have reproduced a simulation from Ref. [9S] using our new collison oper- 
ator, and we find that hyper-viscosity is no longer necessary to obtain steady-state 
fluxes (Fig. 6.7). This can be understood by examining the linear growth rate spec- 
trum of Fig. 6J3 We see that in this system energy diffusion is much more efficient at 
suppressing short wavelength structures than pitch-angle scattering. Consequently, 
no artificial dissipation of short wavelength structures is necessary. 
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Figure 6.8: Linear growth rate spectrum of the entropy mode in a Z-pinch for 
R/L n = 2.0, where R is major radius and L n is density gradient scale length. The 
solid line is the collisionless result, and the two dashed lines represent the result of 
including collisions. The short dashed line corresponds to using only the Lorentz 
operator, while the long dashed line corresponds to using our full model collision 
operator. Both collisional cases were carried out with va = 0.0lv t h,i/R- 
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6.5 Summary 



In Sec. |6.1| we proposed a set of key properties that an ideal dissipation scheme 
for gyrokinetics should satisfy. Namely, the scheme should: limit the scale size of 
structures in phase space in order to guarantee the validity of the gyrokinetic or- 
dering and to provide numerical resolution at reasonable expense; conserve particle 
number, momentum, and energy; and satisfy Boltzmann's if- Theorem. While com- 
monly employed simplified collision operators or hyperviscosity operators may be 
adequate for some calculations [30] , it is important to be able to use the more com- 
plete collision operator described in this paper, which preserves all of these desirable 
dissipation properties. 

In Sec. 6.2 we presented the model collision operator derived in Chapter 5 and 
discussed some of its features that strongly influence our choice of numerical imple- 
mentation. In particular, we noted that local conservation properties are guaranteed 
as long as the (1, vi\, v 2 ) moments of the kp = component of the gyroaveraged 
collision operator vanish. Further, we argued that the collision operator should be 
treated implicitly because in some regions of phase space, its amplitude can be large 
even at very small collisionalities. 

Our numerical implementation of the collision operator was described in Sec. 



6.3| We separate collisional and collisionless physics through the use of Godunov 
dimensional splitting and advance the collision operator in time using a backwards 
Euler scheme. The test particle part of the collision operator is differenced using 
a scheme that possesses discrete versions of the Fundamental Theorem of Calculus 
and integration by parts (upon double application). These properties are necessary 
in order to exactly satisfy the desired conservation properties in the long wavelength 
limit. The field particle response is treated implicitly with little additional computa- 
tional expense by employing repeated application of the Sherman-Morrison formula, 
as detailed in Appendix [F] 

In Sec. 6.4 we presented numerical tests to demonstrate that our implemented 
collision operator possesses the properties required for a good gyrokinetic dissipation 
scheme. In addition to these basic properties, we showed that the implemented 
collision operator allows us to correctly capture physics phenomena ranging from the 
collisionless to the strongly collisional regimes. In particular, we provided examples 
for which we are able to obtain quantitatively correct results for collisionless (Landau 
or Barnes), resistive, and viscous damping. 

In conclusion, we note that resolution of the collisionless (and collisional) 
physics in our simulations was obtained solely with physical collisions; no recourse 
to any form of artificial numerical dissipation was necessary. 
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Chapter 7 

Numerical framework for coupled turbulent transport 
calculations 



7.1 Overview 

As discussed in Chapter 1, any realistic model of turbulent transport and 
heating in hot, magnetized plasmas must account for the interaction of small-scale, 
rapid fluctuations and large-scale, slowly evolving equilibrium profiles. The wide 
range of scales that must be resolved makes direct numerical simulation prohibitively 
expensive. Consider, for instance, the range of time and space scales expected to 
be present in ITER (Tables [TTT~| and 1.2). In order to resolve the turbulent electron 



dynamics, the required grid spacing would be on the order of 10~ 4 centimeters 
perpendicular to the magnetic field in space and 10~ 6 seconds in time. To capture 
the evolution of equilibrium profiles, the simulation domain would have to be on the 
order of 200 centimeters in space and 2 — 4 seconds in time. Assuming several grid 
points are necessary to resolve the smallest scales, this results in 10 6 — 10 7 grid points 
in each of two spatial dimensions (with an additional 10 — 100 along the magnetic 
field) and 10 8 — 10 9 grid points in time. This is in addition to the two dimensions 
of velocity space, each of which requires at least 10 grid points. All told, this comes 
out to approximately 10 25 grid points - a factor of 10 10 larger than the largest fluid 
turbulence simulations possible on today's fastest supercomputers. 

Clearly, the brute-force approach outlined above is not possible on any timescale 
of interest. There is, however, a way forward. The vast separation of time and space 
scales between the fluctuations and the equilibrium can be exploited to allow for 
the separate evolution of the two. The theoretical formalism detailing the process 
of scale separation is described in Chapter [3j While these equations are coupled, 
the time and space scales addressed by each equation are fundamentally different. 
In this chapter, we detail the numerical framework of Trinity [95] . a turbulent 
transport and heating code, which exploits the scale separation of these equations 
to greatly reduce the computational effort necessary to simulate turbulent transport 
and heating. 



The moment equations (3.71) and (3.105) describing the evolution of equilib 



rium density and temperature profiles involve only slowly varying scales in space 
and time. This is by construction - in deriving the equations, we averaged over the 



intermediate time and spatial scales defined by Eqs. (3.48) and 3.19 Consequently. 
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numerical solution of these equations can be achieved using a course space-time grid, 
with spatial variation only in the radial direction. 

In order to evaluate the equilibrium evolution equations, one must specify 
values for the time- and space-averaged turbulent fluxes and heating. Whle the 
averaged quantities evolve on equilibrium scales, the fluxes and heating themselves 
evolve on much faster fluctuation scales. We must therefore conduct turbulence 
simulations on the fluctuation scales and space-time average. Since we only need 
the fluxes and heating on a coarse space-time grid for the equilibrium evolution 
equations, there is no need to simulate the turbulence everywhere in the device, nor 
to simulate it over the entire length of the discharge. Instead, we can simulate the 
turbulence in small regions of space and time (Fig. 1.6) and couple these regions 
together using the equilibrium evolution equations. The only constraint is that the 
space-time domain for the turbulence must be sufficiently large to include the longest 
wavelengths of the turbulence and to reach a steady-state. This clearly provides a 
significant savings in both spatial and temporal resolution. Returning to ITER as 
our example, we see in Table [LT] that the longest turbulent timescale is on the order 
of 10 -4 seconds. Since the equilibrium evolution scale is approximately 1 second, 
we obtain a factor of 10 2 savings by coarse-graining in time. The longest turbulent 
wavelengths perpendicular to the field are on the order of 10 centimeters. Assuming 
several turbulence regions are necessary to sample the device volume, we find little 
savings for a device the size of ITER in the radial direction. There is a way to 
obtain considerable savings in the poloidal direction however, to which we now turn 
our attention. 



7.2 Coupled flux tube approach 

The local, or flux tube, model [11] introduced in Chapter 1 allows for the simu- 
lation of microturbulence in a thin tube, several turbulence decorrelation lengths in 
each dimension, enclosing a single magnetic field line (depicted in Figs. |l.l| and 7.1 ). 



Because the typical spatial scale of fluctuations is much shorter across field lines than 
along field lines, the flux tube is highly elongated along the field line. The critical 
assumption that allows for the reduction of the simulation domain to a single flux 
tube is statistical periodicity. That is, we assume the variation of equilibrium quan- 
tities occurs on such a large scale (relative to the flux tube) that the turbulence is 
homogeneous within the flux tube. We stress that this does not disallow large-scale 
gradients in the problem - only variation of these gradients over the small width 
of a flux tube. The flux tube approximation is thus valid as long as the turbulent 
spatial scale is well separated from all other spatial scales in the problem. This 
means that any effects arising from intermediate spatial scales (such as magnetic 
islands, if present) are not included in the flux tube model. The validity of the flux 
tube approximation in the limit of p* = p/a ~ (k^L)^ 1 <C 1 has been verified nu- 



merically [12]. This is illustrated in Fig. 1.2, where we see that flux tube simulations 



agree very well with global simulations (which do not assume statistical periodicity) 
when the turbulent scale length is much smaller than the equilibrium scale length. 



104 



Figure 7.1: Flux tube from GS2 simulation of the spherical tokamak, MAST. The 
flux tube simulation domain wraps multiple times around the toroidal circumference, 
but covers only a fraction of the anular flux surface it is used to map out (shown in 
light blue). Graphic courtesy of G. Stantchev. 
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In axisymmetric magnetic field configurations, a single flux tube, which com- 
prises only a fraction of a flux surface, can be used to map out the entire flux surface 
(due to the notion of statistical periodicity of the turbulence). This constitutes a 
significant savings in simulation volume, which depends on the toroidal mode num- 
bers of interest in the experiment (in particular, the longest significant wavelength 
present along the field line). We can estimate the savings as follows: in the direction 
perpendicular to the magnetic field line, but contained within the flux surface, the 
flux tube must be large enough to resolve the longest turbulent wavelength. Taking 
k± ~ n^q/r, where is toroidal mode number, q is the safety factor, and r is the 
distance from the magnetic axis to the flux surface of interest, we have 

L± ~ - — ~ ~ — — , (7-1) 

k± n^q n^q 

where Lg is the approximate circumference of the tokamak in the poloidal direction. 
Therefore, the simulation domain of a flux tube covers a fraction of approximately 
l/n^q of a magnetic flux surface. For ITER-like fusion devices, the longest perpen- 
dicular wavelength expected to be important is approximately kj_pi ~ 0.1, corre- 
sponding to a toroidal mode number of ~ 100. Combined with a safety factor 
q > 1, this translates into a savings of factor greater than 100 in simulation volume 
to simulate a single flux surface. 

A single flux surface is not sufficient, however, when evolving radial equilibrium 
profiles. In that flux surface represents a single radial grid point on our 



coarse spatial grid. This is illustrated in Figs. 7.2 and 7.4 In Fig. 7.2, we show a 
poloidal cut of a tokamak and indicate a series of flux surfaces for a representative 
magnetic field configuration. We see that the area of the cross-section accounted for 
by several flux tubes is a small fraction of the total area (as long as the device has 
a sufficiently large radial extent). The actual simulation domain for a single flux 
tube is a rectangular box, which is deformed in physical space as it follows along the 



sheared magnetic field line (Fig. 7.3). Depending on the parallel decorrelation length 



of the turbulence, the flux tube may traverse the toroidal circumference multiple 



times, sampling a poloidal cut of the flux surfae at multiple locations. In Fig. |7.4 
we plot a respresentative radial temperature profile and indicate the coarse grid 
composed of several flux tubes. In each case, notice that the flux tubes have a finite 
radial extent. This is necessary in order to contain several turbulence decorrelation 
lengths. In the limit of an infinitely large device (or, equivalently, infinitely small 
turbulence decorrelation lengths), the radial extent of the flux tubes would shrink to 
a single radial point. There is an intermediate regime in which the flux tubes have 
finite size and, depending on the number of flux tubes and the size of the device, 
may overlap radially. This is not necessarily cause for alarm. While each flux tube 
has finite radial extent, it is representative of a single radial point: equilibrium 
profiles and gradients are taken to be constant across each flux tube. Consequently, 
overlapping flux tubes do not actually sample the same physical space and do not 
lead to double-counting of the turbulence. 
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Figure 7.2: (left): poloidal cross section of a typical tokamak. solid lines indicate 
the shape of magnetic flux surfaces and colored regions indicate a typical portion 
of the tokamak represented by the coupled flux tube approach, (right): cartoon 
illustrating the simulation domain (illustrated in blue) in a poloidal cut at the 
outboard midplane for a single flux tube (representing a radial point, or flux surface) 
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Figure 7.3: Cartoon illustrating the flux tube simulation domain (illustrated in blue) 
for a poloidal cut at the outboard midplane. This flux tube represents the entire 
flux surface, which serves as a radial grid point in our transport equations. 
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Figure 7.4: Cartoon illustrating the portion of the radial temperature profile sampled 
by the use of coupled flux tubes. Each of the blue 'IT shapes represent a flux tube. 
Although each flux tube has finite radial extent, it represents a single radial point 
at the center of its domain. 
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7.3 Normalization of the transport equations 



We start from the gryokinetic transport equations ( |3.71| ) and (3.93) derived 
in Chapter |3j For simplicity, we neglect classical effects and time variation of the 
equilibrium magnetic field (which formally evolves on the equilibrium time scale, but 
in practice evolves on the slower, resistive time scale). As mentioned in Chapter [3j 
we assert that external sources enter our hierarchical equations via the transport 
equations. Consequently, we include general sources in our treatment in this chap- 
ter. We also choose not to calculate h nc , and therefore drop the neoclassical terms 
appearing in our equilibrium evolution equations. However, we do not neglect neo- 
classical effects completely; we include them by using an analytic estimate for the 
neoclassical ion heat flux from Ref. [99], which we add to the turbulent ion heat 
flux appearing in our equations. Since we will not be employing h nc , we henceforth 
drop the t subscript on the turbulent part of h, using h = h t . Also, since we are 
dealing exclusively with equilibrium densities, temperatures, and pressures, we drop 
the nought subscript on these quantities. The resulting equations are: 



dn s 
~dt 

dp s 
dt 



dip d 

dip d 
dVdi) 



dV 

dip 
dV 
dip 



d 3 v (R 2 V(j) ■ Vx) h t 



mv 



R<t>-V X )h 



/rf 3 v^(/? 2 V0-V X ) 

J fM,s 



dF, 



M,s 



dip 



0,s 



(7.2) 
(7.3) 
(7.4) 



J2n s v: u (T u -T s ), (7.5) 



where n s is the equilibrium density, T s is the equilibrium temperature, h s is the 
non-Boltzmann part of the perturbed gyrokinetic distribution function, x is t ne 
generalized potential defined in Eq. (4.5), R is the major radius, <p is the toroidal 
angle, ip = ip p /2ir is a measure of the poloidal flux, 



(j4.5|), R is the major radius, 

is the energy exchange rate 



between species s and species u given in Eq. (3.91 ), V is the flux-tube volume defined 
in Eq. (3.20), and the flux surface average of J 7 , ((J 7 )), is defined in Eq. (??). 



We would like to rewrite the term (i? 2 V0- Vx) in a more enlightening form. To 



begin this process, we need the following identity, derived in Eq. (3.52) by assuming 
an axisymmetric equilibrium magnetic field (B c 



x V0 + RB T V(p): 



v • Vip = -R 2 V<p • (v x B c 



(7.6) 



Using the definition of v x from Eq. (3.47), we have 



v x x B 



R 2 V(p ■ v x x B c 



c|bx V% j x b 

Vx - b ( b • vx 



cR 



B 



T 



RV<P • Vx - ( b • Vx 



(7.7) 
(7.8) 

(7.9) 
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From the gyrokinetic ordering, Vyx <C Vj_X- Consequently, the second term on the 
right-hand side of Eq. (7.9) can be neglected, leaving 



cR 2 V<p • Vx = i? 2 V0 • v E x B 

= -v x -W> 



(7.10) 
(7.11) 



Using the result of Eq. (7.11), and taking advantage of the fact that n s and T s are 



constant on flux surfaces, we have 
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(7.13) 
(7.14) 
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(7.18) 
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(7.20) 



where we have added general external sources, S n and S p , to the equations. 
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The normalized, flux-surface-averaged, nonlinear radial fluxes and turbulent 
collisional heating calculated in Trinity are 



r " s (7 - 21) 



e " s ^W" (7 - 22) 

„» s ///^JLCTiJk^X ( , 23) 



with 

f- = r^V^ (7-24) 
{nTv th ) r pj. 

where the subscript r denotes the reference species in the Trinity calculation, and 
a is a user-specified normalization length. For definiteness, we choose a to be half 
the diameter of the last closed flux surface (LCFS) at the elevation of the magnetic 
axis. We note that p r is defined so that it is a flux surface quantity. Specifically, it 
is given by p r = v thjr /Q a , where Q a = \e\ B a /m r c and B a is the toroidal field on the 
flux surface R a , the average of the minimum and maximum values of R for the flux 
surface of interest. We further define the normalized quantities 

n s = ^~ t = ^ (7.26) 

n 2 

V7 — V7 — V thO,r Pr,0 , / 7 07 -, 

V = av r = Trt (7.27) 

a er 

c a S n a 2 ~ _ a «S P a 2 

O n = 2~ = Tf, 2~ 

VthO,r ^0,r P r fl V t hO,r ^0,r 1 0,r P r n 



VthO,r P r ,0 P 0,s Pr 

where n , r and T 0>r are chosen to be 10 20 m -3 and 1 keV, respectively, and 



_ V t hO,r 2T 0r m r c 3 

= "nT = V ~^7^ -4.57x10 m, (7.30) 

with 5 r chosen to be 1 T. For later convenience, we list here the values of some of 
the normalizing quantities: 



/ II I 

v th0 , r « 4.38 x 10%/-^ m/s (7.31) 



(7-32) 

9.14 V m r y ' 
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where m p is the proton mass and a is the minor radius given in meters. 

Rewriting the transport equations in terms of these normalized quantities gives 
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To proceed, we need some more definitions from Trinity: 
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Using the above result, we get 
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The normalized form of p r is given by 
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B~ a 
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T, 



1/2 



rf + 



8 P Q ' 



(7.40) 



(7.41) 



leading to the following form for the transport equations: 



dn s 



d 



a n dp 



A N 
5? 



+ S n 



(7.42) 




rf + 



<91nT s 
dp 



Q 



N 



(7.43) 



As an aside, we note that in Trinity we choose to keep the density and 
temperature for the reference species equal to unity in each flux tube. This means 
that the physical At and k±p range vary from flux tube to flux tube (since these 
quantities are normalized by v t h, r and p r factors, respectively, which have radial 
dependence). 



7.4 Discretization of the transport equations 

Now that we have a set of normalized equilibrium evolution equations, we pro- 
ceed to discretize them. In doing so, we wish to maximize computational efficiency. 
Primarily, this is achieved by developing an implicit scheme (based on Newton's 
method). While the implicit scheme requires considerably more computational ef- 
fort at each time step than an explicit scheme, it allows for much larger time steps. 
Since the calculation of the steady-state turbulent fluxes and heating at each equi- 
librium time step is by far more expensive than the advancement of the equilibrium, 
the time step size is much more important than the time spent in calculation during 
each step. 

For simplicity, we assume our system consists of electrons and a single ion 
species (which we take to be the reference species), and we use quasineutrality to 
relate the electron and ion densities. Further, we restrict ourselves to the use of a 
three-point spatial stencil. Also, in the interest of notational convenience, we drop 
the tilde on all normalized quantities. 



7.4.1 Particle transport 



We begin by writing the normalized particle transport equation (7.42) for the 
reference (ion) species in the convenient form 



dn 



(|Vpl) 9F 
A dp 



(7.44) 
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where we have defined 



B 2 V/ 2 ' 



(7.45) 



Note that we have switched variables from T to p = nT to reduce the heat transport 



equation (7.43) to an equation involving the evolution of a single variable. A general 



time discretization for Eq. (7.44) takes the form 
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(\Vp\)dF 
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A dp 



+ S n 



, (7.46) 



where At is the transport time step size, and a G [0, 1], with a = corresponding 
to a fully explicit scheme and a = 1 corresponding to a fully implicit scheme. The 
superscripts m and m + l represent the time step. From now on we will drop the 
superscript m wherever it appears; whenever a time superscript is absent from a 
time-dependent quantity, it is understood to be evaluated at time step m. 
Discretizing the spatial derivative using centered differences, we obtain 
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(\Vp\) j F + -F_ 
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(7.47) 

where Ap = the spatial grid spacing, and the subscripts ± indicate 

evaluation at the spatial locations ^±1/2 — ( x j + x j±i)/2, with the subscript j 
denoting the spatial grid index. 



We see that Eq. (7.47) is a nonlinear partial differential equation. We would 



like to treat it implicitly in order to take large transport time steps. This requires 
linearization of the problem. We accomplish this by employing Newton's method, 
in which we expand the m + l time level nonlinear term F about its value at the 
time step m. Keeping terms in this Taylor expansion through linear order, we have 



pm+l ^ 



yoj 



dF± 
dy 



(7.48) 



y=yo 



where y = {{n k }, {pi k }, {p ek }} is a vector containing the values for density and 
electron/ion pressure at each of the spatial grid points, and y is the vector y 
evaluated at time step m. For convenience, we will henceforth drop the y = yo 
specifier on the term dF±/dy. 



Explicitly writing the second term in Eq. (7.48), we have 
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(7.49) 
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Evaluating the partial derivatives of F± in this expression yields 
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(7.50) 
(7.51) 
(7.52) 



efc 



where 5j t k is the Dirac delta function. Substituting Eqs. ( 7.49 )-( 7.52) in Eq. (7.48) 
we obtain the following expression for F™ + 1 
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(7.53) 



The above equation involves derivatives of the particle flux, T, with respect 
to the equilibrium density and pressure at each of the grid locations, n k and p k . 
Unfortunately, this information is not readily available and would be prohibitively 
expensive to compute directly. In order to make calculation of these derivatives 
feasible, we make the assumption that T depends on the {n k } and {p k } only through 
the gradient scale lengths R/L n and Rl p - This assumption is motivated by empirical 
results from both numerical simulation and experiment. 

With this assumption, the derivatives of the particle flux, T, can be written 
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d(R/L p ) ± 

In order to evaluate these expressions, we need discretized forms for R/L n and R/L p , 
as well as estimates for the derivatives of the flux with respect to these equilibrium 
gradients. We defer discussion of the latter issue until later. For the discretization 
of (R/L n ) ± we use 
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(7.56) 



This derivative approximation is accurate to 0[(Ap) 2 ], and the same discretization 
scheme is used for (R/L p ) ± . 

We next compute the discrete derivatives of the equilibrium gradients with 
respect to the equilibrium density and pressure: 
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with a similar expression for the derivative of R/L p with respect to species pressure. 
Using this expression in Eq. (7.53), we arrive at the following: 
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With Eq. (7.58), we can now compute OF /dp: 
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Plugging Eq. (7.59) back into Eq. (7.47), we arrive at the final form for our 
discretized particle transport equation: 
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where 
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7.4.2 Heat transport 

We begin by recasting the normalized heat transport equation 
of species pressure p = nT: 
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To keep notation compact, we define 
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With these definitions, Eq. (7.72) becomes 
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A dp 

Henceforth, we drop the subscript s; whenever a species subscript is not present, 
the subscript s is assumed. 

Our treatment of heat transport in this subsection follows closely our the 
treatment of particle transport in the previous subsection. Our time discretization 
is of the form 
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where a was defined following Eq. (7.46). As before, any time-dependend quantity 
without a time superscript is understood to be evaluated at time step m. 

Linearizing the nonlinear terms via Taylor expansion about density and pres- 
sure at time step m, we again obtain expressions of the form 
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Explicitly writing the second term in Eq. (7.82), we have 
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Evaluating each of the partial derivatives in this expression, we obtain 
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2 Kj A 



dp. 



dp 



d In Tj + dlnp js + 3 d{R/L p ) 3 - 



dp, 



'k 



dp, 



<-'k 



2 k j <9p, 



(7.84) 
(7.85) 
(7.86) 



Substituting Eqs. (7.83)-(7.86) into Eq. (7.82) gives 
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3 djR/Ljj 



(7.87) 
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The expressions for the other nonlinear terms are derived in a similar 
The analogs to Eqs. ( 7.84 )-( 7.86) for H, K, and E are 



dHj 
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dp 



'k 



din Qj 
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dhiQj 



dp 



_ —s 1 djR/L 
2nj 3 (9n^ 
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ik 



dp ek 

dKj 
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djq 

dp 



Hj 



dlnQj _ 1 d(R/L p ) 

. 9p e , Kj dp e . 



K 4 
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dEi 



= K 3 
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dp ek 
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'dlnHj 
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dlnHj 
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2rij _ 



+ 



ik 



2 Pl 



dp 
T , dlnHi 



-6 



<9p 



'efc 
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dp, 



3 5 se m s Z u + d S i7n u Z s 
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- 5 se Z s 3 5 si m s Z u + 5 se m u Z s 
PujZ u — p Sj Z s 2 m s Z u p Uj + m u Z s p Sj 



EjS jk 



EjSjk, 



manner. 

(7.88) 
(7.89) 
(7.90) 
(7.91) 
(7.92) 
(7.93) 
(7.94) 
(7.95) 
(7.96) 
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giving 
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(7.97) 



(7.98) 
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Finally, we consider We have 
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(7.99) 

(7.100) 
(7.101) 
(7.102) 



efc 



which gives us 
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± 



(7.103) 



dp ik ek k dp ek 

At this point in the calculation, we once again make the assumption that T, Q, and 
Ti depend on the {n k } and {p k } only through the gradient scale lengths R/L n and 

With this assumption, the derivatives of the fluxes at grid locations can be 
written 



dT, 



dY j d{RjL n ) 3 



dn k d(R/L n )j dn k 

dTj d(R/L p ) 3 



dp k d(R/L p )j dp k 



(7.104) 
(7.105) 
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with similar expressions for derivatives of the heat flux, Q, and the heating, 7i. We 
have already given discrete forms for R/L n and R/L p at the off-grid locations Xj±i 
in Eq. (7.56). We now give an expression for evaluation at grid locations: 



R 



R f dlnn 

a 



dp 



R rij+i — n 3 -\ 
a 2rijAp 



(7.106) 



This derivative approximation is accurate to 0[(Ap) 2 ], and the same discretization 
scheme is used for R/L p . We next compute the discrete derivative of the equilibrium 
gradients with respect to the equilibrium density and pressure: 
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(7.107) 

with similar expressions for the derivatives of R/L p with respect to species pressure. 



Substituting Eqs (7.106) and (7.107) into Eqs. (7.87), (7.97), (7.98) and (7.99) 



and combining E, G, H, and K results in the following: 
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where 
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(7.109) 
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Next we compute OF /dp: 
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Finally, we add all of these terms up to get 
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where 
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7.4.3 Boundary conditions 

To complete our numerical prescription, we must supply boundary conditions 
at the innermost and outermost radial grid locations. At the outer boundary, which 
corresponds to a location in the fusion device just inside the edge pedestal region, 
we are free to specify density and pressure. The purpose of the simulation is then to 
determine the core temperature as a function of the pedestal density and pressure 
and of the external heat source strength. 

At the internal boundary, which corresponds to the magnetic axis, we use the 
physical constraint that the product of the flux surface area and the flux surface- 
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averaged fluxes is zero: 



lim A | Q t 



(7.129) 



We take the magnetic axis to correspond to the spatial index j = 1/2. Consequently, 



the terms in Eqs. ( 7.46 ) and ( 7.81 ) involving the radial derivative of the fluxes reduce 
to 



(d_F_ 
\dp 



'3/2 



F\ji 



'3/2 



(7.130) 



All other inner boundary terms only involve evaluation at j = 1. However, the 
linearization of the nonlinear terms introduces quantities like d{R/L n )j/dni t . Since 
we have used a three-point, centered stencil for R/L n and R/L p , we would need to 
evaluate the density and pressure inside the inner boundary. To avoid this, we must 
employ an alternate discretization for R/L n and R/L p at j = 1. We choose to use 
a shifted three-point stencil to retain second order accuracy of derivatives: 



R 



R /<91nn 

a 



dp 



R 1 

a rijAp 



-rij + 2n j+1 - -n j+2 . (7.131) 



Computing the derivative of this expression with respect to the density at grid 
points, we have 
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(7.132) 

with similar expressions for pi and p e . 

A consequence of this inner boundary condition is that the coefficients given in 
Eqs. (7.63)-(7.71 ) and (7.116)-(7.128) are modified by taking F_ to be zero. There 



are additional modifications to Eqs. (7.116)-(7.128) due to the new discretization of 
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R/L n and R/L p . The new coefficients are 
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(7.133) 

(7.134) 

(7.135) 

(7.136) 
(7.137) 

(7.138) 

(7.139) 

(7.140) 

(7.141) 
(7.142) 

(7.143) 

(7.144) 

(7.145) 



7.5 Time averaging 

The steady-state turbulent fluxes used in the transport equations are time- 
averaged values. We would like to minimize simulation time by running each set of 
turbulence calculations just long enough to obtain good statistics on the converged 
fluxes. Consequently, we have developed an adaptive time averaging procedure for 
the turbulent fluxes and collisional heating that automatically detects when the 
fluxes have converged. 

For each flux tube simulation, we keep track of the instantaneous time average 
of the fluxes and heating: 



r m = ^r,(At) t , 



(7.146) 
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where m denotes the turbulent time step and At is the size of the turbulent time 
step. Analogous expressions are used for the heat fluxes and heating. While this 
quantity is sufficient for use as the time-averaged flux in our transport calculation, 
we have the additional burden of determining when this flux has converged so that 
we may end the turbulence calculation. In order to accomplish this, we store each 
of the Ti for i = 1 — > m. At each time step we then compute a measure of the rms 
deviation of each of the T* (for % — j — > m — 1, where j is by default m/2, but can 
be specified by the user) from the current value: 



When er is less than a user-specified tolerance, the flux is determined to be con- 
verged, and the turbulence calculation terminates. 

7.6 Quasilinear fluxes 

Since turbulent flux calculations are computationally expensive, we find it 
convenient at times to use quasilinear estimates for the fluxes and heating. We do 
not claim that these estimates are quantitatively correct; we merely use them as a 
computationally inexpensive tool to test our numerical scheme and to gain quick, 
qualitative insight into transport and heating processes. 

To obtain quasilinear estimates for the fluxes and collisional heating, we nor- 
malize the fluxes and heating from linear Trinity simulations by |$| and multiply 
by a factor derived from mixing length theory. The argument goes as follows: one 
expects saturation to occur when nonlinear effects become dominant, i.e. when 
dh/dt ~ v x • Vh. From this balance, we obtain an approximation for the growth 
rate in terms of $: 




(7.147) 




(7.148) 



Defining <3> = (a/p)(e$/T), k = k±p, and 7 = (0/1^)7, we have 



=s>r ; 
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(7.149) 



(7.150) 



with analagous expressions for the heat fluxes and collisional heating. 



7.7 Heat source 
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Taking into account normalizations, we have 
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(7.153) 



For circular flux surfaces, we have 
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We are currently using an analytic source of the form 
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for which 
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(7.156) 



Using this result, we can rearrange Eq. (7.154) to solve for the normalized source 
amplitude A in terms of the input power and the width of the Gaussian power 
deposition profile a: 
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(7.157) 



a£V>(l-exp [-&])' 
where V is given in Mega Watts and rhi is given in units of the proton mass. 



7.8 Trinity simulations 

In this section we present simple tests showing that our implicit transport 
solver is well behaved and that it provides significant computational savings over an 
explicit solver. Furthermore, we present preliminary results from Trinity simula- 
tions using gyrokinetic, turbulent fluxes and heating. These simulations, which are 
the first of their kind, demonstrate that the coupled flux tube approach can rou- 
tinely be used to obtain steady-state equilibrium profiles of density and pressure, as 
well as the corresponding turbulent fluxes and heating. 
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7.8.1 Tests 



The first test we present is intended to demonstrate that the heat transport and 



turbulent heating terms from Eq. (3.105) have been properly implemented and that 
the implicit transport solver algorithm developed in this chapter is well behaved, 
even for multi-channel transport (here we are evolving density and ion and electron 
pressures). We artificially set the temperature equilibration term to zero and use 
the following analytic model for our normalized fluxes and turbulent heating: 

n l/2 

r * = 172 ( 7 - 158 ) 

Pi 

Qs = ^r 2 (7.159) 

Pi 

where all quantities are understood to be the normalized versions defined in Sec. 



7.3 We note that these fluxes do not mesh well with the approximation employed 
in our transport solver algorithm that the fluxes and heating depend primarily on 
R/L n , R/L Ti , and R/L Te . Consequently we are also testing the resiliency of our 
scheme. The resultant transport equations are 

which has the solution p s = F[t + p], with F an arbitrary functional. For our 
initial conditions, we take n{r = 0) = 1 and p s (r = 0) = exp[— p]. Our boundary 
conditions are n(p±, r) = 1 and p s (p±, t) = exp[— (p± + r)], were p± represents the 
the inner and outer radii in the simulation. The solution to this system is n(p, r) = 1 



and p s (p,r) = exp[— (p + r)]. In Fig. 7.5, we show the numerical solution for this 



system, which is in excellent agreement with our analytic prediction. 

Our second test illustrates the superiority of our implicit implentation to an 
explicit scheme when considering fluxes that lead to diffusive behavior. We artifi- 
cially set the turbulent heating and temperature equilibration terms to zero and use 
the following form for the turbulent heat flux (we do not evolve the density in this 
case, so there is no need to define the particle flux): 

Qs = -ir-j ^> (7-163) 

where D is a constant diffusion coefficient. The resultant transport equations are 
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Figure 7.5: Comparison of the analytic and numerical solutions to the system defined 
by Eqs. (7.161) and (7.162) at r = and r = 2. Lines represent analytic solution 



and dots represent numerical solution from Trinity. Here we are showing only the 
ion pressure, but the solution for the electron pressure is identical (and the density 
remains approximately constant in time). Simulation conducted with At = 0.02 
and 16 equally spaced radial grid points (flux tubes). 
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Figure 7.6: Comparison of the analytic solution to the D = 0.1 diffusion equation 
(7.164) at r = (solid line) and r = 2 (dashed line) to the numerical solution from 



Trinity (square dots). Here we are showing only the ion pressure, but the solution 
for the electron pressure is identical. Simulation conducted with At = 0.1 and 16 
equally spaced radial grid points (flux tubes). 



which is simply the diffusion equation. Taking the initial condition of the form 
p s (p,T = 0) = exp[— p 2 /4D] and boundary conditions of the form p s (p±,r 



7.6 



a/1 /r exp \—p\j 4:Dr] , the solution is p s (p, t) = \/l/r exp [— p 2 /ADt] . In Fig. 
we show the numerical solution using Trinity's implicit transport solver. After 
conducting a number of both explicit (a = 0) and implicit (a = 1/2) numerical 
simulations, we find that the implicit scheme gives good results (relative error less 
than 10%) for at least At = 2.0, whereas the explicit scheme is numerically unstable 
for approximately At > 0.02. After taking into account the fact that an additional 
set of flux tube simulations must be run for each transport channel at each trans- 
port time step when running implicitly, we find that the implicit scheme provides a 
savings of a factor of ~ 25 — 50 over the explicit scheme, depending on the number 
of transport channels used (from 1 — 3 currently). 
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7.8.2 Preliminary results 

The simulation results presented in this section are taken from relatively low 
resolution simulations with reduced physics models. They are meant to be demon- 
strations that the coupled flux tube approach detailed in this chapter can be rou- 
tinely used to obtain predictions for steady-state equilibrium profiles and turbulent 
fluxes. All simulations considered here were run with a hyperviscous dissipation 
model [100J that allows us to obtain reasonable, converged turbulent fluxes with 
a relatively coarse spatial grid for the turbulence. For each simulation we used a 
16 x 16 grid in the spatial plane perpendicular to the magnetic field, 26 grid points 
along the magnetic field line, 12 energies, 20 untrapped pitch angles, and a variable 
number of trapped pitch angles (the number of trapped pitch angles in Trinity 
depends on location along the equilibrium magnetic field line. See Chapter 4 for 
more details.). 

We consider a three different systems. All of them have: a single kinetic, 
hydrogenic ion species; electrostatic fluctuations; major radius of 6.2 meters; aspect 
ratio of 3.1; local (Miller) geometry with concentric, circular flux surfaces; fixed 
edge temperature of 4 keV; and external heat input to the ions (via a Gaussian 
deposition profile with a = 0.2/?). The first two systems we consider both have 
adiabatic electrons and 60 MW deposited in the ions from an external source, but 
they have different magnetic field strengths. Evolving only the ion pressure gradient, 
the simulations were run with eight radial grid points and At = 0.02 (At ~ 0.018 
seconds) for 25 time steps. The simulations took approximately 20 minutes each 
on 2048 processors. The steady state ion temperature profiles for the two different 



toroidal magnetic field strengths are shown in Fig. 7.7 As expected, the case with 
the stronger magnetic field leads to higher core temperatures. In Fig. 7J3 we compare 
the B a = 5.3 T result for the ion temperature profile with the same result calculated 
using only neoclassical fluxes (obtained using the analytic expression for the ion heat 
flux from Ref. [99J). We see that in the absence of microturbulence, the core ion 
temperature is well in excess of what is required to ignite a burning plasma. 

The final case we consider has kinetic electrons, a magnetic field strength 
of B a = 5.3 T, and 120 MW external heat source, with 30% going into the ion 
channel. Again using 8 radial grid points, we evolved both the ion and electron 
equilibrium pressure profiles. For our time step, we used At = 0.005 (At ~ 0.004 
seconds) and evolved for 25 time steps. The simulation took approximately one 



hour on 4096 processors. The results, shown in Fig. |7.9| indicate that the use of 
kinetic electrons (instead of the adiabatic electron model used to obtain Fig. 7.7) 
leads to a significant (approximately 65%) reduction in the core ion temperature. 
The fact that the core ion temperature is reduced upon taking into account kinetic 
electron effects is not suprising since the trapped electron mode (TEM) is enabled 
when considering kinetic electrons. However, the large size of the reduction may be 
misleading, since we are employing a relatively coarse grid in phase space. 
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Figure 7.7: Steady-state ion temperature profile for two different values of B a , the 
magnetic field magnitude at the center of the LCFS. As expected, an increase in B a 
leads to an increase in the core temperature. 
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Figure 7.8: Comparison of steady-state ion temperature profiles for simulations 
using turbulent fluxes (solid line) and only neoclassical fluxes (dashed line). Without 
the fluxes arising from microturbulence, core plasma temperatures would easily be 
sufficient to ignite the plasma. 
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Figure 7.9: Steady-state ion and electron temperature profiles for the same system 
used to obtain the B a = 5.3 T plot in Fig. |7.7[ with the exception that here we retain 
kinetic electron effects. Temperature equilibration is strong enough near the edge 
(due to low electron temperature, moderate collisionality, and weak local external 
heating) to keep the ions and electrons in thermal equilibrium, but this is not true as 



we approach the core. Comparing with Fig. 7.7 we see that the core ion temperature 
is significantly decreased by retaining kinetic electron effects. 
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7.9 Summary 



In this chapter, we have detailed a numerical framework for efficiently simu- 
lating turbulent transport and heating in magnetic confinement fusion devices. In 
Sec. 7^2 we introduced the local approximation, which allows for the use of a turbu- 
lence simulation domain consisting of a thin tube encompassing a single magnetic 
field line. Each flux tube is used to map out an entire flux surface, constituting a 
significant saving in simulation volume. These flux surfaces are then used as radial 
grid points in a coarse spatial grid when solving the equilibrium evolution equations 
(3.71) and (3.105) derived in Chapter 3. The steady-state turbulent fluxes and 
heating calculated in each of these flux tube simulations are then time-averaged, 
representing a single step in a coarse, equilibrium timescale grid. 

In Sees. 7.3 and 7.4, we normalized and discretized the equilibrium evolution 
equations. An important consideration in our time discretization was the stiffness 
of the equations, which led us to develop a fully implicit scheme. This was accom- 
plished using Newton's method, in which we expanded the nonlinear terms about 
their values at the previous timestep and kept only terms through linear order. As 
a result, we are forced to evaluate derivatives of the averaged turbulent fluxes and 
heating with respect to the density and pressures at each of the radial grid loca- 
tions. Since this is computationally very expensive, we made the approximation 
that the fluxes depend primarily on gradient scale lengths; the dependence on the 
local density and pressures, as well as the dependence on higher order derivatives, 
is considered to be weak enough so that it can be neglected in taking the flux 
derivatives. 

In Sec. 7.5 we detailed the method by which we obtain numerical time averages 
of the turbulent fluxes and heating. By comparing the running time average to the 
history of accumulated time averages, we defined a criterion that is used to determine 
when the time averaged turbulent fluxes and heating have converged to their steady 
state value. Once they have converged, the flux tube calculation is terminated and 
the time averaged fluxes and heating are passed to the transport solver. 

We described our simple quasilinear flux model in Sec. 7J) and our external 
heating source in Sec. |7.7 Finally, we presented the results of Trinity simulations 



in Sec. 7.8 These results included simple tests showing that the implicit, multi- 
channel transport solver employed in Trinity is well behaved and computationally 
efficient. Additionally, we presented preliminary results from full-volume Trinity 
simulations of the entire discharge of ITER-like plasmas. These simulations, which 
calculated steady state equilibrium profiles and corresponding gyrokinetic, turbulent 
fluxes, constitute the first such simulations ever conducted. Each simulation took 
less than an hour on no more than 4096 processors, making it possible to routinely 
run such simulations in the future. 
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Chapter 8 

Summary and discussion 



In this thesis, we have presented a complete theoretical (Chapter 3) and nu- 
merical (Chapter 7) prescription for studies of the self-consistent interaction between 
turbulence and equilibrium profiles. In order to make such numerical studies fea- 
sible and to ensure that the relevant physics processes are accurately modeled, we 
developed and implemented velocity space resolution diagnostics (Chapter 4) and a 
model physical collision operator for gyrokinetics (Chapters 5 and 6). Combining 
all of these elements, we have developed a new code, Trinity, with which we have 
produced the first ever nonlinear, gyrokinetic simulations of coupled turbulence, 
transport, and heating over a full fusion device volume and discharge time (Chapter 

7)- 

Thus far, the physical systems we have considered have been somewhat sim- 
plified. However, the capability currently exists to do more physically realistic sim- 
ulations, including multiple species, electromagnetic effects, electron scales, and 
general geometry. Consequently, Trinity can immediately be used to explore a 
variety of interesting, experimentally relevant problems. It can be used to conduct 
both qualitative and quantitative studies of possible novel effects of the turbulence- 
equilibrium interaction, such as the formation of internal transport barriers and the 
effect of turbulent heating on the electron-ion temperature ratio. The fact that 
Trinity simulation runtimes are relatively short also allows us to routinely carry 
out parameter scans to study things such as: shaping effects; scaling of performance 
with device size, aspect ratio, and magnetic field strength; and dependence of core 
temperature profiles on edge temperature. 

There are still numerous improvements which could be made to the numer- 
ical algorithms employed by Trinity. Consideration should be given to how one 
can quickly determine whether or not a given nonlinear simulation is stable (below 
the critical gradient threshold) so that runtime is not wasted simulating decaying 
turbulence. Runtime could also be saved by developing a scheme to minimize the 
number of, or utilize use the resources from, idle processors (which appear because 
some flux tube calculations converge faster than others). A preconditioner for the 
profiles calculated using quasilinear or gyrofluid estimates for the fluxes could be 
employed to provide initial profiles to Trinity that will quickly converge to steady 
state. One could explore whether multiple iterations of the Newton solver (instead 
of the single iteration method developed in Chapter 7) allows for the use of larger 



137 



transport time steps and more rapid convergence to steady state. This could also 
lead to the development of an adaptive transport time step. The spatial stencil 
used for finite differences could be widened at essentially no additional computa- 
tional cost. Finally, it may be possible to develop a scheme for evaluating terms 
such as dQ/d(R/L p ) that does not require additional nonlinear flux calculations 
(through the use of quasilinear flux estimates or something similar). 

Possible future directions for improvement to the physics model in Trinity 
include: treatment of the large scale radial electric field, equilibrium flows, and mo- 
mentum transport; inclusion of additional plasma parameter dependencies (such as 
the electron-ion temperature ratio) when approximating the fluxes at the new time 
step in our transport solver algorithm; a more sophisticated (numerical) calculation 
of neoclassical transport and heating effects; treatment of fast particles; more real- 
istic particle, momentum, and heat sources; and treatment of the slow evolution of 
magnetic flux surfaces. Equilibrium shear flows associated with the large scale radial 
electric field profile are believed to play a critical role in the reduction of turbulence 
and formation of the edge pedestal. The development and implementation of mo- 
mentum transport equations would thus allow for quantitative studies of transport 
barrier formation. Inclusion of electron-ion temperature ratio dependence in our 
transport solver would allow us to calculate heating in the homogeneous plasmas 
present in astrophysical systems. In particular, it would allow us to determine the 
ratio of the turbulent energy deposited in ions to the turbulent energy deposited in 
electrons, giving the steady state electron-ion temperature ratio. 

By itself, the first-principles turbulent transport code presented here is not 
sufficient to provide self-contained, comprehensive predictions for the performance 
of fusion devices such as ITER. There are a number of critical physics phenomena 
not currently present in our model, such as equilibrium flows, edge physics, and 
MHD processes. However, it can provide first-principles predictions for core profile 
evolution over a wide range of experimental configurations and plasma parameter 
sets. Furthermore, a code such as Trinity is a necessary component in full-physics, 
predictive simulations of tokamak discharges. Full-physics, predictive simulations 
are a critical component for the fusion program as we develop ITER and look beyond 
to the next generation of fusion devices. 
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Appendix A 
Geometry 

A.l General geometry 

Our development closely follows that of Ref . [Uj . Since the divergence of the 
magnetic field is zero, one may use a Clebsch formulation [101] : 

B = Vax Vt/>. (A.l) 

To represent an equilibrium magnetic field composed of closed surfaces, it is sufficient 
to define |lUlj : a = <fi — q(ip)9 — u(i/j, 9, <p) and ip — \ff. Here, 9 and are 
the physical poloidal and toroidal angles, respectively, \I/ = (2n)^ 2 f v dr~B ■ V9 is the 
poloidal flux, g(^) = d^ T /d^, ^ T = (2tt)~ 2 f y drB ■ V0 is the toroidal flux, and dr 
is the volume element. The quantity v should be periodic in 9 and 4>. 

It is convenient to define a new angle ( = <fi — v. With these definitions, 
Eq. A.l becomes 

B = W x V(g0-C), 

where the subscript on B is included to emphasize that we are concerned with the 
equilibrium, unperturbed magnetic field. The field lines are straight in the ((, 6) 
plane, and are labeled by a. Useful coordinates are therefore (p,a,9), where p(^) 
determines the flux surface, a chooses a field line in that surface, and 9 measures 
the distance along that field line. 

In an axisymmetric system, one may also represent the magnetic field as 

B = /(^)V0 + x V0, (A.2) 

where = RBt- We will find it useful to take advantage of this representation, 
although not necessary. 

In the ballooning or field-line following limit, we assume that the perturbed 
quantities vary as 

A = A(9) exp (iS) 

where b ■ VS = 0. This takes into account the fact that the perturbations tend to 
be slowly varying along the field line, and allows for rapid variation across the field 
line 0. 

The latter condition implies 

(Va x W) • VS = 

which, in turn, implies S = S(a, To make contact with the ballooning approxi- 
mation and with field-line following coordinates, one may choose S = no (a + q9o), 
where uq is some (large) integer, and #o is the familiar ballooning parameter which, 
in field-line- following coordinates, determines k x through the relation k x = —kgs9o- 
Here, s = p/q(dq/dp), and p is an arbitrary flux surface label. 
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A. 2 Operators and arguments 



In general, we wish to simulate the nonlinear electromagnetic gyrokinetic equa- 
tion in the ballooning, or field-line-following, limit. We choose a field-line-following 
representation [llj, which has the advantage that the nonlinear terms are easy to 
evaluate and are independent of the details of the magnetic geometry. Further de- 
tails may found in Ref. Below, we focus on the linear terms, which may be 
affected by the geometry. 

Effects of the magnetic geometry in this limit enter through only a small num- 
ber of terms, regardless of whether one proceeds with a moment-based approach [TT] , 
a 5f approach |102(, I103j . or a gyrokinetic approach |1041 [5l]. Consider, for example, 
Eqs. (23-24) of Ref. [8]: 



g = h 



B dji 



, , v ± \VS\\ ( 2 v\\ 7 \ ,v ± \VS\ fv ± \VS\ 



, (A.3) 



and 



— i ( to — uJd + iv ||b ■ V ] h 



exp (-iL)st f 



8F B x VS • VF 



de 



\VS\v ± j fv ±\VS\ \ 

n J 



'-qijj ) + qa— — —Ji 
c / c 



(A.4) 



Here, cu d = VS ■ B x ( mv 2 b • Vb + fiVB + gV$ J /{mB Q). The notation is 

explained in Ref. [8]. Note that the unperturbed magnetic field B = B (9). 

These equations, together with Maxwell's equations, describe the linear prop- 
erties of a wide range of microinstabilities. In the limit of large toroidal mode number 
n , only the following components of these equations depend on 9: b • V, (VS"! 2 , 

B x ^b • Vbj ■ VS, (B x VB ) ■ VS, and B (9). To perform volume integrations 
and flux surface averages in the nonlinear simulations, it is also necessary to have 
the Jacobian J and |Vp| as functions of 9. We now consider the terms individually. 

To make our normalizations clear, we treat the uj* term in detail. The oj* term 
may be written as 



B x VS- VF A . c „ 
~ % B mn qX = - m °B- X 



b x V (a + q9 ) ■ VF ( 



where 



This, in turn, is 



c / c 



c A 
Bo 



b x V (a + q9 ) ■ VF 



-m -x (b-V«x V*j — = -m cx— , 



140 



where we have assumed that F = F (\l/). 

We now introduce normalizing quantities. Lengths are normalized to a, which 
we choose to be half the diameter of the last closed flux surface (LCFS), measured at 
the elevation of the magnetic axis. The magnetic field is normalized to the toroidal 
field on the flux surface at R a , (B a = I(^)/R a ) where R a is the average of the 



minimum and maximum of R on the flux surface and I(ip) is as used in Eq. (A. 2). 
Time is normalized to a/v t , where v t = ^T/rrii. Thus, for example, V = (l/a)VAr 
and \1/ = a 2 B a ^/ N . Perturbed quantities are scaled up by a/p ia , where p ia = v t /fl a 
and Q a = \e\B a / \rriic). The perturbed field is normalized by Ti/\e\, so that, for 
example, xn = (\ e \x/Ti)(a/ pi a )- [Here, we consider only the one-species problem. 
The generalization to multiple species is straightforward.] Finally, we introduce an 
arbitrary flux surface label p, normalized so that p = at the magnetic axis and 
p = 1 at the LCFS. Note that the Larmor radius pi should not be confused with the 
flux surface label p. Upon adopting these normalizations, one finds 

„<9F n cT p ia A dF dp p ia v t dF 0/s 
OW a 2 eB a a Op dw N a 2 op 

which serves to define kg = (n /a)dp/d^/ ^. In the high aspect ratio, zero /3, circular 
flux surface limit, kg = n^q/r. For the case in which there is a background density 
gradient, one finds 

p ia V t dF A F p ia V t .,. . „ a (PiaVt\ 

-ikgpia — - — — xn = i{k e p ia )XNTi—, = i{k d p ia )XNF — — =- 

a 2 Op {L n ) N a 2 L n V a 2 / 

in which the dimensionless quantity (L n )~^ = —(l/n)dn/dp, and may also be writ- 
ten as L n /a. With the specified normalizations for time, space, and perturbed 
quantities, the factor pi a Vt/a 2 scales out of the gyrokinetic equation. Compare, for 



example, the cu* term with the first term in Eq. (A. 4) 



• , , ,'PiaVt 

lUJri = lUNflN 



a 2 



The factor in parentheses is common to all terms in the equation, and does not 
appear in any other form. It may therefore be considered to be arbitrary. 

In the uj* term, note that ke is multiplied by pi a , confirming that it is natural 
to consider perpendicular gradients normalized by the gyroradius pi a rather than to 
the minor radius a, as expected in the ballooning or field-line-following limit. 



To summarize, upon adopting the above normalizations, the term in Eq. (A.4) 
in field-line-following coordinates becomes 

.BoXVS-VFo A . — „ (ptaV t \ ., ldFo. „ (Pia,Vt\ /a r- \ 

1 — B^n — qx = iuj * nXnFo {-*-) = - lkep -Y ^p- XNFo {-*-) (A ' 5) 

Note that a;*jy = —kePiai^-/ Fo)(dF /dp) is dimensionless, independent of 9, and 
related to the dimensional u;* by cu* = u>*NVt/a- 

We now turn to the b • V operator. We begin by using the B field in the form 



of Eq. flA.lD to find a: 

da 

B • V0 = W x W ■ V0— 

do 
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which implies 



a 



d6- 



B - V0 



(A.6) 



For an axisymmetric B field, this integral may be evaluated with the use of Eq. (A. 2 ). 
In this case, the b • V operator may be explicitly evaluated. It is 



b ■ Vh(8) 
which serves to define 



B ■ W9 dh 

Bn 86 



1 N 



aB 



N 



da 
86 



IV 



,dh 
86' 



b- V 



N 



B 



N 



da 
dd 



-l 



IV 



iVS 



N 



(A.7) 
a/qRo, where 



In the high aspect ratio, zero j3, circular flux surface limit, \h • V 

Ro is the major radius at the center of the flux surface. 

Next, we consider the V B part of the u d operator. This term is given by 



h 



2 VLBl 



B o xV5 -V.S 



PiaVt 



-LN 



N 



2 dtf 



N 



B% dp 



b x V N B N ■ Vat (a + q6 ) 



The module released here produces the factors in square brackets, i. e. 

2 d$ N 



2 d^ N ~ {0) 
wvb = ^2"— i— b x V N B N ■ V N a and u^' B 

JJjy Up 



B 2 N dp 



bx V N B N -V N q. (A.8) 



In the high aspect ratio, zero /3, circular flux surface limit, uj^ b — 2a /Ro (cos 6 + s6 sin I 



and uj- 



(0) 



-2 (a/R )s sin9. 



The curvature drift is nearly the same as the V-B drift, except that v]_ — > 2vi. 
and the fact that there is an additional component of the curvature drift given by 



, 47T/l 

qbJ 



bxVp-V5 



PiaVt 



kePi 



h N vl N 



1 dtf 



N 



B% dp 



b x V N (3 a -V N (a + qd 

(A.9) 



The module released here produces the factors in square brackets, i.e., 



1 dV N ~ 

CJ K = CJVB + T^T— 1— b X v NPa ■ VatO, 



a;- 



(0) 



(A.10) 



Here, /3 a = 8irp/Bl. Note that a perpendicular gradient of (3 a gets no contribution 
from the gradient of the magnetic field, since B a is a constant. 

We do not explicitly consider the remaining component of Ud, proportional to 
V$o- To the extent that the electrostatic potential is constant on a flux surface, 
this term may be specified using the information provided by the module. 

To summarize, in field-line-following coordinates, the term in Eq. (A.4) that 
is proportional to uJd is given by 



iujdh 



he pi, 



h 



N 



PiaVt 



u VB + 6 J® B ) + vjL (u K + 6 4 0) ) 
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The form of Eqs. (A.3- A~4~| ) and of the gyrokinetic Maxwell's equations [8] 
(not shown) guarantees that \VS\ always appears as the square, IVS"] 2 . In general 
geometry, this term may be written 



IVS1 



n| 
a 2 



Vat (a + q9 



N 



dp 



\{V N a + 9 V N q) ■ (V N a + 9 V N q)\ 



The module released here produces the factors (gi, #2, #3), where 

2 



|V5| 2 = k% \g x + 29 g 2 + 9 2 g 3 \ = k 2 



N 



dp 



I Vjvtt • V^a + 29 V N a ■ V N q + 9^V N q ■ V N q\ 



In the high aspect ratio, zero /3, circular flux surface limit, gi 



(A.ll) 

1 + s 2 9\ g 2 = -9s, 



and gs = s 2 . Note that IVS 1 ] 2 typically appears with a factor of 1/fi 2 , which is not 
included in Eq. (A.ll). 

The remaining quantities are straightforward. The variation of the unper- 
turbed magnetic field along the field line is reported by the module as B^, with 



B N {9) = B (9)/B a . 



(A.12) 



The quantity |Vjvp| is also reported by the module, and is unity in the high aspect 
ratio, zero (3, circular flux surface limit. 

For numerical applications, it is sometimes necessary to choose the field line 
coordinate so that (b • V)jv is constant. This choice allows the straightforward eval- 
uation of terms proportional to \k\\ \ in the transform space. Thus, we use (p,a,9') 
coordinates, where 9' is the equal arc periodic coordinate defined by 



9'(9) = 2ttL n (9)/L n (7c) 



7T 



(A.13) 



-i 



between —it and ir, and Ln{9) = J_ 7T d9 • Vj . In this coordinate system, the 
coefficient of the parallel gradient operator of Eq. (A. 7) becomes 

b- v)' = 2n/L N (n). 



(A. 14) 



The Jacobian is Jn = (d^x/dp) (L n /2ttB n ) . With these definitions, the flux sur- 
face average of a quantity T is defined to be 



/ r Jjv d9' da 



J Jtv dO' da 

The normalized area of the flux surface is A^v = 27r( |Vtvp| ) / JdO' . 



The fi eld-line variation of the quantities cuy b and cjyg [Eq. (A. 8)], uj r a nd uj^ ' , 



,(o) 



[Eq. l^lCf], ( 9u92 , 93 ) [Eq. B N m [Eq. £l2|], and . Vj^ [Eq. ^14f], 

together with the quantities | vnp\, dp/dty n and d(3/dp are the outputs of this geom- 
etry module. These coefficients contain all of the geometric information necessary 
for numerical calculations of high-n microstability and turbulence in axisymmetric 
toroidal configurations with nested magnetic surfaces. 
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A. 3 Module details 



Input numerical equilibria may be specified in numerous ways, as documented 
in the module. Interfaces to direct and inverse Grad-Shafranov equilibrium solvers 
are available. These include using output from several equilibrium codes in use 
in the fusion community, such as TOQ [1051 M, EFIT [TUT], VMOMS [T08] . 
JSOLVER jTUn], and CHEASE [110], as well as the local equilibrium model of 
Ref. UH]. 

Here, we describe our implementation of the Miller local equilibrium model [lllj 
for completeness. This model extends the usual zero-beta, high-aspect ratio equi- 
librium to arbitrary aspect ratio, cross section and beta, and allows one to consider 
geometric effects on microinstabilities in a controlled way. 

The shape of the reference flux surface and its perpendicular derivative are 
specified in the (R, Z) plane by 

R N (6) = R 0N (p)+p cos [9 + 5(p) sine], (A.15) 

Z N (0) = K(p)p sin 6. (A. 16) 

Here, R N = R/a, etc., R n{p) = RoN{Pf) + R'oN d P, S(p) = 5(p f ) + 5' dp, n(p) = 
K i.Pf) + ^ dp, and pj denotes the flux surface of interest. In the remainder of this 
section, the "N" subscripts will be dropped, since no ambiguities will arise. 

As emphasized in Ref. [lllj . the actual shape of neighboring flux surfaces 
(p Pf) is n °t determined by Eqs. (A.15) and (A.16). Instead, this is determined 
by solving the Grad-Shafranov equation in the neighborhood of pf. As noted by 
Mercier and Luc |112j . one may find this solution provided R{6), Z{6), B p {6), p'(pf), 
and I'(pf). One additional piece of information is required to determine either the 
safety factor q or / dp. Finally, the normalization of the magnetic field determines 

In our implementation, we take q to be an input parameter, and upon noting 
that § add = —lisq, use it to define d^> /dp from Eq. (A.6): 

dV 



d6 



(we x vp • v^) 



-i 



(A.17) 



dp 2nq J R? 

[For numerical equilibria, d^ /dp may be calculated directly, and this expression 
defines the safety factor.] The poloidal magnetic field B p (pj) is specified by 

dV \Vp\ 



B n 



dp R 



where |Vp| may be found from Eqs. (A.15) and (A.16). 

The remaining steps may be used with the Miller local equilibrium model or 
with full numerical equilibria. We allow arbitrary values of dp/ dp and s by using the 
mercier expressions to find VS [1111 11121 1113] . As noted in Ref. the result 

is exactly equivalent to the generalized s — a analysis of Greene and Chance |114] . 
To proceed, we define 



MP) 



de 



R 2 



+ 



B P R 2 



B(e) 



de 



(B P RY 
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C{9) 



de 



sinw + R/R c 
B V R" 



V0 x w • 

where u(9) is the angle between the horizontal and the tangent to the magnetic 
surface in the poloidal plane, and R c is the local radius of curvature of the surface 
in the poloidal plane. If we define A — § ■ ■ ■ , etc., it can be shown that 

q dp 2irq dp v ' 

where the primes denote derivatives with respect to \I> '. Thus, one may specify any 
two of p', I f , and s. This freedom is a direct consequence of the two free functions 
in the Grad-Shafranov equation. 

It can also be shown |113[ llllj that the perpendicular component of the gra- 
dient of a is given by 

V«-e* = |W| (AI' + Bp + 2C) . 



The parallel component of the gradient of a may be easily found from Eq. (A. 6). 
With Va in hand, the remainder of the calculation is straightforward. We note that 
V-B may also be calculated using the mercier formulas; our treatment is the same 
as can be found in Refs. [11 lj and |113j . To wit, the perpendicular component is 

^ „ * Br, ( B v . I 2 sin u \ 
VB ■ e* = + Rp' - — , 



B \R c R 3 B p j 

and since B(pf) does not depend on p' or the component of VB along the field 
line depends on neither quantity. 

The expressions for s and the gradients of a and B make it clear that once 
the safety factor, the shape of the flux surface, and B p are determined (either from 
a numerical equilibrium or from the local equilibrium), one may vary p' and s in- 
dependently to find a family of solutions, all of which satisfy the Grad-Shafranov 
equation. This flexibility allows one to carry out the Greene-Chance kind of analysis 
for microinstabilities. Such an analysis simplifies the interpretation of the numerical 
calculations, since all other parameters can easily be held fixed. 

Within the context of the local equilibrium model |lllj . one may also vary 
individual shape parameters one at a time, to explore the dependences in a controlled 
fashion. 

The eleven dimensionless parameters that determine the local MHD equilib- 
rium in this implementation of the Miller model are summarized in Table I. 
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In addition to these eleven parameters, there are two normalizing dimensional 
parameters, a and B a . In all, there are two more parameters than are found in 
Ref . [11 lj . We include the additional parameters to allow a somewhat more natural 
correspondence between reported equilibria and the input variables. We emphasize 
that there is nothing "extra" in our implementation of the model as result; it is 
only modestly easier to use for some applications. For example, our choice of the 
normalization of the magnetic field (B a ) is the vacuum magnetic field at -R geo , the 
center of the LCFS. This quantity is the most commonly reported magnetic field 
value. By allowing Rq to be specified separately, we also make it conceptually easier 
to separate the effects of Shafranov shift from the derivative of the Shafranov shift. 
The inclusion of the normalized minor radius as a separate variable is a natural 
choice as soon as one allows for separate specification of R and -R geo - 

The starred quantities (pf, s, and d(3/dp) may be specified when reading in 
numerical equilibria. Values of the latter two quantities that are different from the 



actual equilibrium values are incorporated by using Eq. (A. 18) to define I'. 

Finally, when using numerically generated equilibria, the module allows one 
to choose from the most common definitions of p, such as the normalized poloidal 
or toroidal fluxes, the horizontal minor radius, etc. The user may also provide his 
or her own definition of p by supplying a simple function. 
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Appendix B 

Landau damping of the ion acoustic wave 

We consider the collisionless ion acoustic wave in slab geometry with adiabatic 
electrons. The gyrokinetic equation for this system has the particularly simple form 

dt dz T dt K J 

Changing variables from h to g = {f\) and assuming solutions of the form 

g = g(v) e i k u z - iuJt ) , (B.2) 

we obtain 

( u -kv)g = kv^^F M , (B.3) 

where I am using v = v\\ and k = k\\ for convenience. Neglecting FLR effects and 
assuming quasineutrality gives 

(w -kv)g = kvr^- J d 3 v'g(v'). (B.4) 

Defining 

PCX) 

g(v) = 2tt / v ± dv ± g(v) (B.5) 
and integrating over the perpendicular velocities in the gyrokinetic equation yields 

(u-kv)g(v)=kF(v), (B.6) 

where 

F(v) = vr^^e'^, (B.7) 
V2irv t 



J dv'g(v')- (B.8) 



Following the analysis of Case and van Kampen, we see that this equation has 
solutions of the form 



V— l — + \(k,u)5(u-v) 
u — V 



B.9) 



g(v) = F(v) 

with m = f , provided that A is chosen to satisfy the condition 

n 1 = J dv'g(v') = V J dv'^l + \(k,u)F(u). (B.10) 
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A general solution is given in the form 

/OO POO 
/ C(k,u)g k>u (v)e ik ^ ut Uk du, 
-OO J — oo 

where C(k, u) is determined by the initial condtion 



(b.ii; 



f{z,v,0) 




C(k,u)g ku (v)e tkz dk du. 



Taking the inverse Fourier transform of the above expression gives 



F(k,v) 



C{k,u)g k Jv)du, 



where 



f(k,v) = — I f(z,v,0)e-* kz dz. 

Z7T 



Plugging the expression (B.9) for g into the initial condition (B.13) yields 

C(k,u) 



T{k, v) 
F(v) 



V 



u — v 



-du + \(k,v)C(v). 



(B.12) 



(B.13) 



(B.14) 



(B.15) 



We now have two equations, (B.10) and (B.15), for two unknowns (A and C) 



In order to solve this linear system, it is convenient to define some new notation. 
Any square integrable function H can be written 



H(q) 



K(p)e ipq dp. 



(B.16) 



We define the positive and negative frequency parts of H as 

/■±oc 

H ± (q) = ± / K(p)e^dp, 
Jo 



(B.17) 



so that H = H + + H_. Further we define the function = H + — H . It can be 
shown that if* has the alternate form 



HJv)=V- 



1 r° H(v') 



7T2 ./_«, V'-V 



dv' . 



With these definitions in hand, we rewrite eqns (B.10) and (B.15) as 



F(v) 



-iriF^u) + \F(u), 

(A + ttz) C+{v) + (A - tti) G-{v). 



(B.18) 

(B.19) 
(B.20) 



Eliminating A gives an expression involving C + and C_: 

F(k, u) = (n 1 + 2niF + (u))C + (k,u) + (n 1 -2mF-(u))C-(k,u). (B.21) 
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The transform T can also be broken down into negative and positive frequency parts 
to give two separate equations. 



F±{k,u) = (n 1 ±27iiF±(u))C±(u) 
These can then be used to construct C(k,u): 



C{k,u) 



+ 



T-{k^u) 



ni + 2itiFj r {u) rii — 2niF(u 



(B.22) 



(B.23) 



Substituting the expressions (p.9|) and (p.23|) for g and C into the equation 

T-(k, u) 



(B.ll ) for / gives 




+ 



rii + 2iiiF + {u) ni — 2-7iiF-(u) 



F(v) 



V- 



u — v 



+ \{k, u)5{u — v) 



Mz-ut) dk du _ 



We can use the identity 



F±{k,u) 



1 

27 



- ikz 'dz' 



to rewrite eqn (B.24) in the more convenient form 

5 + {u — v') 



f(z,v,t) 



V- 



5-(u-v') 

ni + 2iTiF + (u) n\ — 2-niF^[u 
1 



2tt 



F(v) 



u — v 



+ \(k, u)5{u — v) 



iH*-z'-ut) dz f dv f dk du _ 



Now we pick an initial condition of the form 

f(z,v,0) = f(v,0)e ik ° z , 

which gives 



f{z,v,t) = e ife °(*-**) (m + mF^v)) 



/+M) 



+ 



/-M) 



rii + 2iriF + (v) ni — 2niF_{v 



+ V 



F(v) 



U(u,0) 



+ 



/>,0) 



u — v \ nx + 2iriF + (u) ni — 2niF_(u 



D ik (z-ut) du 



(B.24) 



5±(u-v')f(z',v',0)dv' (B.25) 



(B.26) 



(B.27) 



(B.28) 
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Appendix C 

Proof of the iif-Theorem for our model collision operator 

In the case of the expansion / = Fq + 5f about a Maxwellian the entropy 
generation by like particle collisions takes the form 



dS 
~dt 



d 
dt 



flnfdv dr 



fC[fF ] dvdr>0, 



(C.l) 



where we use the compact notation / = 5f /Fq. The statement of the if-theorem is 



that the right-hand side of Eq. (C.l ) is nonnegative and that it is exactly zero when 
Sf is a perturbed Maxwellian. 

We represent / as a Cartesian tensor expansion (or equialently spherical har- 
monic expansion) in velocity space: 



f(r,v) = f (r,v)+v- fi(r,v) + R[f](r,v) 



(C.2) 



where R[f] comprises the higher order terms. It is then possible to recast the 
statement of the if-Theorem in terms of this expansion using linearity of the model 
collision operator C [Eq. (5.13)], orthogonality of the expansion and the fact that 
spherical harmonics are eigenfunctions of the Lorentz operator L. By construction, 
R[f] satisfies J R[f]F dv = and J vR[f]F dv = 0, from which it follows that R[f] 
does not contribute to the field-particle part s of the model operator: U \R\f ] Fn] — 
and Q[R[f)F ) = 0. Substituting Eq. joi} into the right-hand side of Eq. (fcjj, 



where the operator C is given by Eq. (5.13), and integrating by parts those terms 



involving derivatives of R[f], we find that they all give nonnegative contributions, 
so we have 



J fC[fF ]dv>a + a l , 



where 



00 



0"! 



f C[f F ]dv, 

v ■ frCiv ■ fiF \dv. 



(C.3) 

(C.4) 
(C.5) 



In order to prove the if-theorem, it is now sufficient to show that o~ > and a% > 0. 



Starting with ctq and using Eq. (5.13), we integrate over angles and use the 



identity v^veFq = —(d/dv)(v 5 i , ^F ) to express the term containing Q: 



-27T / /, 



d_ 
1 dv 



v\F ^-[f -^Q[f F 



dv 



dv. 



(C.6) 
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Using the aforementioned identity again in the expression for Q [Eq. (5.12)] and 
integrating by parts where opportune, we get 



dfo 
dv 



v vwFndv 



5 771 I 

v i/\\rodv 



dv 



v^ueFqcLv. 



(C.7) 



It is easy to see from the Cauchy-Schwarz inequality that 



(C. 



Using this in the second term of Eq. (]C.7|), we infer 



dfo 



dv 



v v\\F a dv 



\ ~ J v\F Q dv I J v 6 v E F dv^J = 0, 



(C.9) 



where to prove that the right-hand side vanishes, we again used the identity v a veFq = 
— (d/dv)(v 5 v\\F ) and integrated by parts. Thus, we have proved that er > 0. 

Turning now to 0\ [Eq. flC.5 )], using Eq. (5.13), and integrating by parts where 
opportune, we get 



(Ji= / (v ■ f!) 2 u D F dv 



1 



3v. 



th 



2.1 \dv 

2 







xx ■ fxu s Fodv 



v- fx) v u\\F dv 



x 2 u s Fodv, 



(CIO) 



where we have used the standard notation that x = v/v t h and x = v/v^. Integrating 
over angles and using the simple identity a- f vvdQ = (47r/3)a, where v = v/v and 
a is an arbitrary vector, we have 



0"! 



^(f\fi\ 2 x'u D F dx+ 1 - 



d_ 
dx 



xfi 



x 4 u\\Fodx 



fix A v s F dx 



(C.ll) 



X^VeFndx 



Once again applying the Cauchy-Schwarz inequality, we find that 



f x x A v s F Q dx 



< 



fi 



x A is s F dx / x 4 is s F dx. 



(C.12) 
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Using this in the last term in Eq. (C.ll), we get 



47T 



'tli 



fl 



x AvF^dx 

2 



d_ 
dx 



xfi 



(C.13) 



x A is\\F dx 



where Av = vd — v s . Upon using the identity 2x 3 Az/F = {d / dx){x i v\\Fo) in the 
first term of the above expression and integrating the resulting expression by parts, 
we finally obtain 





S/i 




<9x 



x & v\\F dx > 0. 



(C.14) 



We now consider whe n these inequalities becomes equalities, i.e., when the 
right-hand side of Eq. (C.l) is zero. Firstly, this requires dR[f]/d£ = and thus 
R[f] = 0, so / contains no 2nd or higher-order spherical harmonics. Secondly, o"o = 
if either f is independen t of y or we have equality in the invocation of the Cauchy- 
Schwarz inequality [Eq. jC.8 )], which occurs if f oc v 2 . Similarly a% = iff fi is 
independent of v. Thus, the right-hand side of Eq. (C.l) vanishes iff / oc l,v,v 2 , 



l.C 



fFo is a perturbed Maxwellian. 
This completes the proof of the if-theorem for our model operator. 
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Appendix D 

Gyroaveraging our model collision operator 



To transform the derivatives in Eq. (5.13) from the original phase-space co- 
ordinates (r, v,£, i?) to the new coordinates (R, v , £, we require the following 
formulae: 



3_ 
d_ 



d_ 
dv 
d_ 

di 
d_ 



R 



+ 



R 



+ 



R 





(- 


-p 

V 


{or 


t 


„. I 






( d 




\dR 



d_ 
OR 



(D.l) 
(D.2) 
(D.3) 



where p = b x v±/Q. In Fourier-transformed perpendicular guiding center variables, 
we can replace in the above formulae (d/dR) v — > ik, where k = k±. It is convenient 
to align (without loss of generality) the $ = axis with A;, so we have v± ■ k = 
k±v yl — £ 2 cos-# and p ■ k = —k\_vyl — £ 2 sint?. Using the above formulae, we 



gyroaverage the Lorentz operator in Eq. (5.13): 

1 d 



(L[h k \) 



2d£ 



dhk 



v 2 (i+e 

An 2 



-kj_h k , 



(D.4) 



where we have used {pp) R = (v±v±) R = (1/2)1. Note that both the terms con- 
taining £ and t? derivatives in the original operator Eq. (5.9) produce non-zero gy- 



rodiffusive contributions [the second term in Eq. (D.4)]. Another such gyrodiffusive 
term, equal to —v\\ [v 2 {l — £ 2 )/4£T] k\ h k , arises from the energy-diffusion part of the 
test-particle operator in Eq. (5.13). Collecting these terms together and defining 
the thermal Larmor radius p = Vth/Q, we arrive at the gyrodiffusion term in Eq. 

It remains to gyroaverage the field-particle terms. For the energy-conservation 



(5.19) 



term [Eq. (5.12)] we have 

(e ik - p Q[e~ ik - p h k }) = (e* k - p ) Q[e~ ik - p h k ] 



(D.5) 



where 



Q[e- lk - p h k ] = 

J v 2 v E (e~ ik - p ) h k dv I j v 2 (v/v th ) 2 u E F dv. 



(D.6) 



Note that the d integration in Q only affected e p , hence the above expression. 

Using the standard Bessel function identity [TT7] J e MSlIl1? dd = 2nJ (a), we find 

o 
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(e lfc ' p ) = Jo(a), where a = k±v±/fl. Substituting this into Eqs. (D.5) and (D.6), 
we arrive at the energy-conservation term in Eq. (5.19), where expression in the 



right-hand side of Eq. (D.6) is denoted Q[h k ] [Eq. (5.22) . 

The momentum-conserving terms are handled in an analogous way: details 
can be found in Appendix B of Ref. [IS], where a simpler model operator was 
gyroaveraged. 
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Appendix E 

Comparison of our collision operator with previous model operators 

In order to compare and contrast with previously suggested operators that do 
include energy diffusion, we first rewrite in our notation the operator derived by 
Catto and Tsang[25] [their Eqs. (14) and (16)], 

Cct [6f] = u D L[6f] + ~ (\^ n F Q ^- S -f\ + ^ v ■ I vu s 6fdv 

v 2 dv\2 dvF J n % J 

2Fhf^_S\ r v 2 ' 

3n \v 2 h 2 / •/ < ih 



+ — (T-nl / —VE^fdv. 



This operator, whilst it conserves particle number, momentum and energy, neither 
obeys the if- Theorem nor vanishes on a perturbed Maxwellian. 

The latter point can be demonstrated most easily by letting Sf = x 2 F Q , where 
x = v/vth- This 5f is proportional to a perturbed Maxwellian with non-zero 5n and 
ST. We can then evaluate the test-particle and field-particle parts of the operator 
to find 

1 d fl 4 _ d Sf\ Id 



v 2 dv 



and 



I 



4 °° l~2 
x 2 UESfdv = —j^z \ x 6 UEe~ x dx = \l — uqv. (E.3) 

V 71 " J V 7T 



Substituting into Eq. (E.l), we get 



CcT[Sf] = -x 2 u E F + ^ (x 2 - uF , (E.4) 



which is non-zero despite 5f being a perturbed Maxwellian. 



In order to show that the ii-theorem can be violated by the operator Eq. (E.l ), 



let us consider a perturbed distribution function of the form 5f = x 3 F . Then 



CcT[Sf] 



(E.5) 



so the entropy generation is, 

= " // % Cct [6f] dvdr = ~^ (32 _ 21v/2_) vV < °' (E - 6) 



where V is the volume of the system. The above expression is negative, which breaks 
the if -theorem and produces unphysical plasma cooling for the particular form of 
the perturbed distribution function that we have examined. 



155 



The second case we examine here is the sequence of operators derived by 
Hirshman and Sigmar.|24j The general operator proposed by these authors is given 
in their Eq. (25). In our notation, we rewrite here the N = 0, I = 0, 1 restriction of 
the like-particle form of their operator with Au set to for simplicity (this does not 
affect the discussion that follows): 



C HS [Sf}=u D L[Sf] 



1_ d_ 

v 2 dv 



1 4 P 9 




QM 



J (E.7) 



2vU 

+ vd — n — Fq, 



where U and Q are defined by Eqs. (5.11) and (5.12). The primary concern here 



comes from the angle averaging operation in the energy- diffusion part of the opera- 
tor. Firstly, the energy diffusion only acts on the spherically symmteric (in velocity 
space) part of the perturbed distribution function. However, there is no reason why 
there cannot arise perturbations that have very large energy derivatives but angle- 
average to zero (for example, 8f oc £). Clearly, such perturbations will not damped 
correctly. Secondly, upon conversion to gyrokinetic coordinates and gyroaveraging 



see Sec. 5.3 and Appendix D), the operator becomes 



Cffi^GK^fc] = Vd{v) 

J (a) d 



1 9 dh k 



(l-£ 2 ) 
2 d£ y q J <9£ 



+ 



v 2 dv 



i 



4 

J {a)h k 



2u D 



-1 

v±Ji(a)U± [h k ] +v\\J (a)U\\ [hi 



-F + ve^t J {a)Q[h k ]F , 



'th 



th 



(E.8) 



where the conservation functionals U±, U\\ and Q are the same as defined in Eqs. 
( 5.20 )-( 5.22). The immediatly obvious problem is that the angle averaging has 
introduced two new Bessel functions into the energy diffusion term. The energy 
diffusion is therefore supressed it by one power of k±p in the limit k±p ^> 1, while 
it is precisely in this limit that we expect the small-scale structure in the velocity 
space to be particularly important. [T51 [T4"] This means that the energy cutoff in 
phase space is artificially pushed to smaller scales and one might encounter all the 
problems associated with insufficient energy diffusion .[22] 

While, for the reasons outlined above, we expect the Hirshman-Sigmar oper- 
ator not to be a suitable model for collisions, we would like to note that for many 
purposes the Hirshman-Sigmar operators are superior to the model operator we pre- 
Taken as a sequence, they provide a rigorous way of obtaining 



sented in Sec. 5.2 



classical and neoclassical transport coefficients to any desired degree of accuracy, 
and it is relatively easy to solve the Spitzer problem for them, while the Spitzer 
functions for our operator are hard to find analytically. 
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Appendix F 

Sherman-Morrison formulation 

The repeated application of the Sherman-Morrison formula considered here is 
an extension of the scheme presented in Tatsuno and Dorland. [97] Throughout this 
calculation, we adopt general notation applicable to both Eqns. (6.23) and (6.24) 
and provide specific variable definitions in Table F.l Both Eqns. (6.23) and (6.24) 
can be written in the form 

Ax = b. (F.l) 

Because (U) and (£) are integral operators, we can write them as tensor products 
so that 

A = A Q + u ® v + ui ® Vi + u 2 ® v 2 . (F.2) 

We now define 

Ax = A Q + u ® v , A 2 = A 1 + u 1 ®v 1 , (F.3) 

so that 

(A 2 + u 2 <g> v 2 ) x = b. (F.4) 

Applying the Sherman-Morrison formula to this equation, we find 

v 2 • y 2 



Y2 



-Z2, 



1 + v 2 • z 2 
where A 2 y 2 = b and A 2 z 2 = u 2 . 

Applying the Sherman-Morrison formula to each of these equations gives 

vi • yi 



(F.5) 



y% = yi - r— — : — wi (F.6) 

1 + Vi • Wi 

z 2 = zi - — wi, (F.7) 

1 + Vi • Wi 

where Aiyi = b, AiWi = Ui, and AiZ\ = u 2 . A final application of Sherman- 
Morrison to these three equations yields 

v • Yo 



yi = yo - 

Wi = Wo 

Zi = z - 



1 + v • s 
v • w 



1 + v • s 
v • z 

-s 



so 



so 



(F.8) 
(F.9) 
(F.10) 



1 + v • s 

where A y = b, A s = u , A w = ui, and A z = u 2 . 

We can simplify our expressions by noting that v 0j i,2 and u 0j i,2 have definite 
parity in v\\. A number of inner products then vanish by symmetry, leaving the 
general expressions 



y2 = yo - 
z 2 = z - 



vp • yo 

1 + v ■ s 
1 + v • s 



s - 
so - 



vi ■ yo 

1 + Vi • w 

Vl • Zp 
1 + Vi • w 



w 
w . 



(F.ll) 
(F.12) 
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variable 


L 


D 


A 


1 - At (L + U L ) 


1 - At (D + U D + E) 


X 


h** 


h n+l 


b 


h* 


h** 


A 


1 - AtL 


1 — AtD 


v 


U D Vj_Ji 


—Auv±Ji 


Vi 


VdV\\J 


— Avv\\ Jo 


v 2 


v\\ (electrons) 
(ions) 


u E v 2 J 


u 


-At\ /d u 


-At\ /d u 


Ui 


-Atvx/du 


-At\i/d u 


u 2 


—AtVpV\\/d q (electrons) 
(ions) 


-Atv 2 /d q 


d u 


/ d 3 v v D v^F G 


f d 3 v AvvjF 


d q 


<e/2 


j d 3 v u E v 4 F 



Table F.l: Sherman-Morrison variable definitions for Lorentz and energy diffusion 
operator equations 
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Appendix G 

Compact differencing the test-particle operator 

In this appendix, we derive a second order accurate compact difference scheme 
for pitch-angle scattering and energy diffusion on an unequally spaced grid. The 
higher order of accuracy of this scheme is desirable, but it does not possess dis- 
crete versions of the Fundamental Theorem of Calculus and integration by parts 
when used with Gauss-Legendre quadrature (or any other integration scheme with 
grid spacings unequal to integration weights). Consequently, one should utilize this 
scheme only if integration weights and grid spacings are equal, or if exact satisfaction 
of conservation properties is not considered important. 



For convenience, we begin by noting that Eqns. (6.23) and (6.24) can both be 
written in the general form 



H (Gti)' + S = H (G'ti + Gh") + S, 



(G.l) 



where: for the Lorentz operator equation (6.23) we identify H — 1, G — (1 — £ ) /2, 
S = Ul[H\ — k 2 v 2 vd (1 + £ 2 ) H/AQq, and the prime denotes differentiation with 



respect to and for the energy diffusion operator equation (6.24), we identify 
H = l/2v 2 F , G = v\\v 4 F , S = U D [h] + E[h\ - k 2 v\ (l-f)V^o» and the 
prime denotes differentiation with respect to v. Here, the h we are using is actually 
normalized by Fq. 

Employing Taylor Series expansions of h, we obtain the expressions 



and 



til 



ti 



cL (h 



(h i+ i - hi) + 51 (hi - hi 



hi 



- h (hi 



6-: 



hi. 



+ 0[S 2 



5 + 5_ (5+ + 5.) 



-h? + 0[S 2 



(G.2) 



(G.3) 



where i denotes evaluation at the velocity space gridpoint x«, and S± = |xj±i — x,| 
(here x is a dummy variable representing either £ or v). In order for the h'( expression 
to be second order accurate, we must obtain a first order accurate expression for h'" 



in terms of hi, h\, and h". We accomplish this by differentiating Eqn. (G.l) with 
respect to x: 

'8h'\ 

■ = H' (Gti)' + H (Gti)" + S' (G.4) 

c 



dt 



and rearranging terms to find 



til 



dhj 
dt 



c 



HiG, 

- Hi (G'lti + 2GX 



Hi (G% + Gih'l) 

s: +o[5], 



(G.5) 
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where, unless denoted otherwise, all quantities are taken at the n + 1 time level. 
Plugging this result into Eqn. G.3 and grouping terms, we have 



S- -5, 



h' 



u 1 r^ 1 u r^ 11 

^7 _ Mi^i - 



<L (h i+ i - hi) - 5+ (hi - hi-i) 



At 



S' 



(G.6) 



0[5 2 



5+5- (5 + + 5J) 

where //j = 1 + (5+ — 5-) (H[Gj + 2ifjG-) /ZHGi, and we have taken (dh/dt) c 
(h n+1 — h n )/At. Using Eqn. (G.4) and the above result in Eqn. (G.l), we find 



dh 

di 



o 



hi ( Hid 



5+ + 5- 

3/jLi 



^7 - H i^i ~ tiiGi 



+ 
+ 



HiGi / <L (h i+1 - hi) - 5+ (hi - hi-i) 



5+ + 5. 
3HiGi 



5+5^ (5+ + 5-) 



(G.7) 



At 



+ S' 



+ S t + 0[5 2 



This is the general compact differenced form to be used when solving Eqns. (6.23) 
and fltT24ft . 

In order to illustrate how compact differencing affects the implicit solution 
using Sherman-Morrison, we present the result of using the particular form of S for 
energy diffusion in Eqn. (G.7): 

hr 1 



h? 



hi+\ 



At 



5+ (5+ + 5_) 



ZHiGi 
Hi 



hi-i 



2 HiGi 



5- (5+ + 8 J) v fM 
hi I 2HiGi 



5+G 



5-5, 
At 



[Li 



+ (5+ - 5-) Q - 5+5.Ku t 



(G.8) 



h n 
n i+l 



5- 



5+ (5+ + <L 
5+ -8 



h n - 
n i-l 



5+5- 



5- (5+ + 5-) 
+ U\\V\\ + U+V± + U q q + 0[5 2 



where 



<Ti = (5+ - 5-) /3/ii 

k = k\i (i - e) m 

Q = HiG'i - Ui (l/At - H[G'i - HiG- + Kv s ) 
v s = v\\v 2 th /2v 2 
A = A + a A' 

U±,\\,q = u 0,l,2 



V\ 



V 0,l,2 
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with u and v given in Table F.l 



We see that the only significant effects of compact differencing on numerical 
implementation are: modification of h** in Eqn. (6.24) to reflect the h n terms on the 
right-hand side of Eqn. (G.8); and modification of the U\\, U±, and U q terms that 



appear in Sherman-Morrison 
term. 



(u 0i i,2 from Appendix A) to include an additional all' 
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Appendix H 

Sample Trinity input file 

! case with 2 species, nonlinear fluxes, 4 keV edge temperature 
&nt_params 



! geometry 
rad_out = 0.8 
rmaj = 6.2 
aspr =3.1 
bmag = 5.3 



outer rad. bound (normalized by minor radius, a) 
major radius, R (in meters) 
aspect ratio (R/a) 

B-field amplitude at center of LCFS (in Tesla) 



! species parameters 

ntspec = 2 ! number of species to evolve in transport equations 
qi = 1 ! ion charge (in units of proton charge) 

mi = 1 ! ion mass (in units of proton mass) 



! time advancement 
ntstep = 25 
ntdelt = 0.005 
subcycle = .false, 
nsub = 1 
impfac =0.5 



number of transport time steps 
transport time step size 

set to T to subcycle temperature equilibration 
number of subcycles per time step 
time-centering (0=explicit, l=implicit) 



! fluxes 
grad_option 



'tgrads' 



ql_flag = .false. 
model_flux = .false. 
include_neo = .true, 
dfprim = 1.0 

dtprim = 1.0 
pflx_min =1.0e-5 

qflx_min = 5.0e-4 



assume fluxes depend on ion and electron 
pressure grads 

T for quasilinear estimate for fluxes 
T for offset linear estimate for fluxes 
T to include neoclassical ion heat flux 
amount by which to perturb dens grad when 
obtaining estimate for flux derivative 
amount by which to perturb temp grad 
minimum particle flux to use if GS2 
calculation does not converge 
minimum heat flux 



temp_equil = .true. 



! F to neglect temperature equilibration 
! between species 



! initialization 

rln =1.0 ! initial R/Ln for density profile 

rlti =5.0 ! initial R/LTi for ion temperature profile 
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rite =5.0 
nedge = 1.0 
tiedge =4.0 
teedge =4.0 

! sources 
densin =0.0 
powerin = 120.0 
src_ratio =0.5 
psig = 0.2 
nsig =0.2 

/ 



initial R/LTe for electron temperature profile 
fixed dens at outer rad. bound (in 10"{20}/m"{3}) 
fixed ion temp at outer rad. bound (in keV) 
fixed electron temp at outer rad. bound (in keV) 



amplitude of particle source 

heat source power (in MW) 

fraction of powerin going into ions 

width of Gaussian heat source profile 

width of Gaussian particle source profile 
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